Regression Analyses
Drivers of Maladaptation
an1 <- lm(mal ~ host, data=n1)
summary(an1)
##
## Call:
## lm(formula = mal ~ host, data = n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53932 -0.27745 -0.03932 0.26068 0.72255
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.27745 0.06669 4.160 7.2e-05 ***
## hostC 0.26187 0.08040 3.257 0.00158 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3592 on 91 degrees of freedom
## (53 observations deleted due to missingness)
## Multiple R-squared: 0.1044, Adjusted R-squared: 0.09457
## F-statistic: 10.61 on 1 and 91 DF, p-value: 0.001582
Effects of Maladaptation
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)
# main effect models
an1.1 <- glm(no_tim ~ c.mal + host + ln.vol_cub + con_tim + cn2,
family=quasipoisson,
data=n1, subset=no_tim > 4)
summary(an1.1)
##
## Call:
## glm(formula = no_tim ~ c.mal + host + ln.vol_cub + con_tim +
## cn2, family = quasipoisson, data = n1, subset = no_tim >
## 4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.49224 -0.83543 -0.09667 0.62544 2.36662
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.207594 1.213446 -0.171 0.8652
## c.mal -0.028035 0.298378 -0.094 0.9257
## host1 -0.045736 0.358165 -0.128 0.8992
## ln.vol_cub 0.167263 0.073955 2.262 0.0304 *
## con_tim 0.051004 0.028218 1.807 0.0798 .
## cn2 0.001773 0.035684 0.050 0.9607
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.90865)
##
## Null deviance: 35.583 on 38 degrees of freedom
## Residual deviance: 29.055 on 33 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
an1.2 <- glm(no_tim ~ c.mal + host + ln.vol_cub + con_tim,
family=quasipoisson,
data=n1, subset=no_tim > 4)
summary(an1.2)
##
## Call:
## glm(formula = no_tim ~ c.mal + host + ln.vol_cub + con_tim, family = quasipoisson,
## data = n1, subset = no_tim > 4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.51373 -0.80168 -0.06667 0.59580 2.36304
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.149383 0.952609 -0.157 0.8763
## c.mal -0.005711 0.286871 -0.020 0.9842
## host1 -0.042412 0.154675 -0.274 0.7855
## ln.vol_cub 0.164847 0.071695 2.299 0.0276 *
## con_tim 0.056262 0.024962 2.254 0.0306 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.8610003)
##
## Null deviance: 36.470 on 39 degrees of freedom
## Residual deviance: 29.299 on 35 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
an1.3 <- glm(no_tim ~ host + ln.vol_cub + con_tim,
family=quasipoisson,
data=n1, subset=no_tim > 4)
summary(an1.3)
##
## Call:
## glm(formula = no_tim ~ host + ln.vol_cub + con_tim, family = quasipoisson,
## data = n1, subset = no_tim > 4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.51720 -0.80017 -0.06795 0.59200 2.36087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.14857 0.93844 -0.158 0.8751
## host1 -0.04130 0.14220 -0.290 0.7731
## ln.vol_cub 0.16473 0.07046 2.338 0.0251 *
## con_tim 0.05629 0.02458 2.290 0.0280 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.8370613)
##
## Null deviance: 36.47 on 39 degrees of freedom
## Residual deviance: 29.30 on 36 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
# final model
an1.f <- glm(no_tim ~ ln.vol_cub + con_tim,
family=quasipoisson,
data=n1, subset=no_tim > 4)
summary(an1.f)
##
## Call:
## glm(formula = no_tim ~ ln.vol_cub + con_tim, family = quasipoisson,
## data = n1, subset = no_tim > 4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.51357 -0.79264 -0.04083 0.53831 2.37497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.23088 0.88708 -0.260 0.7961
## ln.vol_cub 0.17084 0.06667 2.563 0.0146 *
## con_tim 0.05505 0.02397 2.297 0.0274 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.8169915)
##
## Null deviance: 36.470 on 39 degrees of freedom
## Residual deviance: 29.371 on 37 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
##### With Neg-Freq Dep Mal Calc -- no difference
an1.1 <- glm(no_tim ~ con_tim + ln.vol_cub,
family=quasipoisson,
data=n1, subset=no_tim > 4)
summary(an1.1)
##
## Call:
## glm(formula = no_tim ~ con_tim + ln.vol_cub, family = quasipoisson,
## data = n1, subset = no_tim > 4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.51357 -0.79264 -0.04083 0.53831 2.37497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.23088 0.88708 -0.260 0.7961
## con_tim 0.05505 0.02397 2.297 0.0274 *
## ln.vol_cub 0.17084 0.06667 2.563 0.0146 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.8169915)
##
## Null deviance: 36.470 on 39 degrees of freedom
## Residual deviance: 29.371 on 37 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)
#### for new maladaptation model -- from p = 0.03 to 0.056
an2.1 <- glm(no_art_5 ~ c.nf_mal_k0.5 + ln.vol_cub + c.con_art5*host,
family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.1)
##
## Call:
## glm(formula = no_art_5 ~ c.nf_mal_k0.5 + ln.vol_cub + c.con_art5 *
## host, family = quasipoisson, data = n1, subset = no_tim >
## 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5517 -1.6891 -0.7854 0.7949 8.0273
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.39575 1.07406 -3.162 0.00216 **
## c.nf_mal_k0.5 0.40763 0.23405 1.742 0.08511 .
## ln.vol_cub 0.46044 0.08406 5.478 4.13e-07 ***
## c.con_art5 0.06330 0.02624 2.413 0.01793 *
## host1 0.41085 0.21784 1.886 0.06262 .
## c.con_art5:host1 -0.05911 0.03704 -1.596 0.11416
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.568432)
##
## Null deviance: 652.94 on 92 degrees of freedom
## Residual deviance: 414.74 on 87 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
####
an2.1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5*host,
family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.1)
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5 * host,
## family = quasipoisson, data = n1, subset = no_tim > 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6751 -1.6447 -0.7855 0.8684 7.8261
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.07216 1.02646 -2.993 0.0036 **
## c.mal 0.49991 0.23689 2.110 0.0377 *
## ln.vol_cub 0.43164 0.08042 5.367 6.55e-07 ***
## c.con_art5 0.06731 0.02611 2.578 0.0116 *
## host1 0.47820 0.21732 2.200 0.0304 *
## c.con_art5:host1 -0.06740 0.03701 -1.821 0.0720 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.346001)
##
## Null deviance: 652.94 on 92 degrees of freedom
## Residual deviance: 407.37 on 87 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an2.2 <- glm(no_art_5 ~ c.mal + ln.vol_cub + host * c.con_art5,
family=quasipoisson, data=n1, subset=no_tim>1)
summary(an2.2)
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + host * c.con_art5,
## family = quasipoisson, data = n1, subset = no_tim > 1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.8730 -1.9463 -0.5069 0.9554 7.5987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.83738 1.27432 -1.442 0.15454
## c.mal 0.44057 0.33068 1.332 0.18780
## ln.vol_cub 0.33915 0.09909 3.423 0.00112 **
## host1 0.19255 0.35264 0.546 0.58707
## c.con_art5 0.05499 0.02912 1.888 0.06382 .
## host1:c.con_art5 -0.03657 0.05358 -0.683 0.49755
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 6.056629)
##
## Null deviance: 431.96 on 65 degrees of freedom
## Residual deviance: 324.94 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an2.3 <- glm(no_art_5 ~ mal + host + ln.vol_cub + host*c.con_art5 + cn_rat,
family=quasipoisson, data=n1)
summary(an2.3)
##
## Call:
## glm(formula = no_art_5 ~ mal + host + ln.vol_cub + host * c.con_art5 +
## cn_rat, family = quasipoisson, data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.8894 -1.5379 -0.5226 0.9420 7.5316
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.69695 1.49892 -1.799 0.07580 .
## mal 0.57739 0.25635 2.252 0.02708 *
## host1 0.76222 0.40843 1.866 0.06572 .
## ln.vol_cub 0.45973 0.08489 5.416 6.39e-07 ***
## c.con_art5 0.08046 0.02808 2.865 0.00534 **
## cn_rat -0.04820 0.04701 -1.025 0.30829
## host1:c.con_art5 -0.06674 0.04068 -1.641 0.10484
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.285246)
##
## Null deviance: 624.76 on 85 degrees of freedom
## Residual deviance: 361.47 on 79 degrees of freedom
## (60 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an2.4 <- glm(no_art_5 ~ mal * host + ln.vol_cub, family=quasipoisson, data=n1)
summary(an2.4)
##
## Call:
## glm(formula = no_art_5 ~ mal * host + ln.vol_cub, family = quasipoisson,
## data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.1205 -1.7340 -0.8982 0.8801 7.5285
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.13525 1.00302 -3.126 0.0024 **
## mal 0.51650 0.31353 1.647 0.1031
## host1 0.64591 0.29033 2.225 0.0287 *
## ln.vol_cub 0.41083 0.07754 5.298 8.54e-07 ***
## mal:host1 -0.29777 0.49148 -0.606 0.5462
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.511588)
##
## Null deviance: 652.94 on 92 degrees of freedom
## Residual deviance: 437.29 on 88 degrees of freedom
## (53 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
# main effect models
an2.1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + cn2 + c.con_art5*host,
family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.1)
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + cn2 + c.con_art5 *
## host, family = quasipoisson, data = n1, subset = no_tim >
## 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.7302 -1.5775 -0.5168 0.9412 7.4701
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.28015 1.48490 -1.536 0.12864
## c.mal 0.56452 0.25661 2.200 0.03074 *
## ln.vol_cub 0.45829 0.08452 5.422 6.22e-07 ***
## cn2 -0.05462 0.04758 -1.148 0.25450
## c.con_art5 0.07966 0.02794 2.851 0.00556 **
## host1 0.80262 0.40694 1.972 0.05207 .
## c.con_art5:host1 -0.06513 0.04052 -1.607 0.11196
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.261526)
##
## Null deviance: 624.76 on 85 degrees of freedom
## Residual deviance: 360.09 on 79 degrees of freedom
## (8 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an2.f <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5*host,
family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.f)
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5 * host,
## family = quasipoisson, data = n1, subset = no_tim > 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6751 -1.6447 -0.7855 0.8684 7.8261
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.07216 1.02646 -2.993 0.0036 **
## c.mal 0.49991 0.23689 2.110 0.0377 *
## ln.vol_cub 0.43164 0.08042 5.367 6.55e-07 ***
## c.con_art5 0.06731 0.02611 2.578 0.0116 *
## host1 0.47820 0.21732 2.200 0.0304 *
## c.con_art5:host1 -0.06740 0.03701 -1.821 0.0720 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.346001)
##
## Null deviance: 652.94 on 92 degrees of freedom
## Residual deviance: 407.37 on 87 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
# partial r-sq
an2.1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1, subset = host=="C")
summary(an2.1) # 37.6% R-sq
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson,
## data = n1, subset = host == "C")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6270 -1.5586 -0.7073 0.7373 7.9524
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.05242 1.26538 -2.412 0.0189 *
## c.mal 0.71362 0.34010 2.098 0.0401 *
## ln.vol_cub 0.42820 0.09911 4.320 5.95e-05 ***
## c.con_art5 0.07210 0.02852 2.528 0.0141 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 6.078791)
##
## Null deviance: 512.65 on 63 degrees of freedom
## Residual deviance: 311.59 on 60 degrees of freedom
## (9 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
n1.p <- n1[complete.cases(n1$no_art_5, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5), ]
n1.p$no_art_5 <- sqrt(n1.p$no_art_5)
# to drop levels properly, keep vars of interest
n1.p <- n1.p[,c("no_art_5", "c.mal", "ln.vol_cub", "c.con_art5", "host", "no_art_5")]
an.full <- lm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
data=n1.p, subset = host=="C")
an.art_mal <- lm(no_art_5 ~ ln.vol_cub + c.con_art5,
data=n1.p, subset=host=="C")
an.art_vol <- lm(no_art_5 ~ c.mal + c.con_art5,
data=n1.p, subset=host=="C")
an.art_con <- lm(no_art_5 ~ c.mal + ln.vol_cub,
data=n1.p, subset=host=="C")
an.mal <- lm(c.mal ~ ln.vol_cub + c.con_art5,
data=n1.p, subset=host=="C")
an.vol <- lm(ln.vol_cub ~ c.mal + c.con_art5,
data=n1.p, subset=host=="C")
an.con <- lm(c.con_art5 ~ c.mal + ln.vol_cub,
data=n1.p, subset=host=="C")
r.art_mal <- resid(an.art_mal)
r.art_vol <- resid(an.art_vol)
r.art_con <- resid(an.art_con)
r.mal <- resid(an.mal)
r.vol <- resid(an.vol)
r.con <- resid(an.con)
an.p_mal <- lm(r.art_mal ~ r.mal)
summary(an.p_mal) # r-sq = 5.91 %
##
## Call:
## lm(formula = r.art_mal ~ r.mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9844 -0.7116 -0.0408 0.5926 3.6861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.776e-17 1.436e-01 0.000 1.0000
## r.mal 9.371e-01 4.210e-01 2.226 0.0297 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.149 on 62 degrees of freedom
## Multiple R-squared: 0.074, Adjusted R-squared: 0.05906
## F-statistic: 4.955 on 1 and 62 DF, p-value: 0.02967
an.p_vol <- lm(r.art_vol ~ r.vol)
summary(an.p_vol) # r-sq = 33.7 %
##
## Call:
## lm(formula = r.art_vol ~ r.vol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9844 -0.7116 -0.0408 0.5926 3.6861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.327e-17 1.436e-01 0.000 1
## r.vol 5.473e-01 9.519e-02 5.749 2.95e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.149 on 62 degrees of freedom
## Multiple R-squared: 0.3477, Adjusted R-squared: 0.3372
## F-statistic: 33.05 on 1 and 62 DF, p-value: 2.945e-07
an.p_con <- lm(r.art_con ~ r.con)
summary(an.p_con) # r-sq = 7.11 %
##
## Call:
## lm(formula = r.art_con ~ r.con)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9844 -0.7116 -0.0408 0.5926 3.6861
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.174e-18 1.436e-01 0.000 1.0000
## r.con 1.143e-01 4.739e-02 2.413 0.0188 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.149 on 62 degrees of freedom
## Multiple R-squared: 0.08584, Adjusted R-squared: 0.07109
## F-statistic: 5.822 on 1 and 62 DF, p-value: 0.0188
# herbivore abundance
anh <- glm(no_herb_5 ~ host*c.con_art5 + cn2 + ln.vol_cub,
family=quasipoisson, data=n1)
summary(anh)
##
## Call:
## glm(formula = no_herb_5 ~ host * c.con_art5 + cn2 + ln.vol_cub,
## family = quasipoisson, data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0557 -1.7191 -0.4770 0.5317 6.5428
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.301610 0.919498 -4.678 7.11e-06 ***
## host1 0.499453 0.276353 1.807 0.07301 .
## c.con_art5 0.050709 0.019107 2.654 0.00894 **
## cn2 -0.036775 0.028114 -1.308 0.19315
## ln.vol_cub 0.609810 0.057103 10.679 < 2e-16 ***
## host1:c.con_art5 -0.008761 0.026901 -0.326 0.74520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 4.263334)
##
## Null deviance: 1314.29 on 136 degrees of freedom
## Residual deviance: 518.29 on 131 degrees of freedom
## (9 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
# on arthropod diversity (5 mm)
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)
## for neg freq dep mal
an3 <- glm(rich_art_5 ~ c.nf_mal_k2.0 * host + ln.vol_cub,
family=quasipoisson, data=n1, subset=no_tim > 0)
summary(an3)
##
## Call:
## glm(formula = rich_art_5 ~ c.nf_mal_k2.0 * host + ln.vol_cub,
## family = quasipoisson, data = n1, subset = no_tim > 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5840 -0.9553 -0.3483 0.7788 2.3267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.72255 0.70002 -3.889 0.000195 ***
## c.nf_mal_k2.0 0.45476 0.19918 2.283 0.024827 *
## host1 -0.06453 0.16319 -0.395 0.693484
## ln.vol_cub 0.35269 0.05493 6.421 6.67e-09 ***
## c.nf_mal_k2.0:host1 -0.64648 0.45049 -1.435 0.154820
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.524261)
##
## Null deviance: 234.49 on 92 degrees of freedom
## Residual deviance: 140.23 on 88 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
##
an3 <- glm(rich_art_5 ~ c.mal*host + ln.vol_cub + host * c.con_art5,
family=quasipoisson, data=n1, subset=no_tim > 0)
summary(an3)
##
## Call:
## glm(formula = rich_art_5 ~ c.mal * host + ln.vol_cub + host *
## c.con_art5, family = quasipoisson, data = n1, subset = no_tim >
## 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3865 -1.0153 -0.2121 0.7193 2.0200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.75709 0.68817 -4.006 0.000131 ***
## c.mal 0.72859 0.21402 3.404 0.001010 **
## host1 -0.02279 0.17607 -0.129 0.897303
## ln.vol_cub 0.35623 0.05413 6.581 3.51e-09 ***
## c.con_art5 0.05687 0.01870 3.041 0.003126 **
## c.mal:host1 -0.98839 0.39786 -2.484 0.014921 *
## host1:c.con_art5 -0.08287 0.02862 -2.896 0.004791 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
##
## Null deviance: 234.49 on 92 degrees of freedom
## Residual deviance: 122.51 on 86 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an3 <- glm(rich_art_5 ~ host*c.mal + host+ ln.vol_cub + host*c.con_art5,
family=quasipoisson, data=n1, subset=no_tim_gs > 0)
summary(an3)
##
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + host + ln.vol_cub +
## host * c.con_art5, family = quasipoisson, data = n1, subset = no_tim_gs >
## 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3865 -1.0153 -0.2121 0.7193 2.0200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.75709 0.68817 -4.006 0.000131 ***
## host1 -0.02279 0.17607 -0.129 0.897303
## c.mal 0.72859 0.21402 3.404 0.001010 **
## ln.vol_cub 0.35623 0.05413 6.581 3.51e-09 ***
## c.con_art5 0.05687 0.01870 3.041 0.003126 **
## host1:c.mal -0.98839 0.39786 -2.484 0.014921 *
## host1:c.con_art5 -0.08287 0.02862 -2.896 0.004791 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
##
## Null deviance: 234.49 on 92 degrees of freedom
## Residual deviance: 122.51 on 86 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
contrasts(n1$host) <- c(0,1)
an3.1 <- glm(rich_art_5 ~ host*c.mal + host+ ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1, subset=no_tim>1)
summary(an3.1)
##
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + host + ln.vol_cub +
## c.con_art5, family = quasipoisson, data = n1, subset = no_tim >
## 1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5165 -1.1008 -0.1576 0.9107 2.0497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.50369 0.82816 -3.023 0.00368 **
## host1 0.61763 0.25225 2.448 0.01729 *
## c.mal -0.83967 0.61093 -1.374 0.17442
## ln.vol_cub 0.28991 0.06540 4.433 4.02e-05 ***
## c.con_art5 0.04197 0.01804 2.326 0.02341 *
## host1:c.mal 1.54806 0.67284 2.301 0.02490 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.446096)
##
## Null deviance: 153.997 on 65 degrees of freedom
## Residual deviance: 90.206 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an3.2 <- glm(rich_art_5 ~ host*c.mal + host+ ln.vol_cub + host*c.con_art5 +
cn_rat,
family=quasipoisson, data=n1, subset=no_tim_gs > 3)
summary(an3.2)
##
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + host + ln.vol_cub +
## host * c.con_art5 + cn_rat, family = quasipoisson, data = n1,
## subset = no_tim_gs > 3)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6129 -0.9568 -0.2084 0.6939 2.2196
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.655505 2.317628 -0.283 0.7791
## host1 -0.390998 0.604090 -0.647 0.5221
## c.mal -0.093110 1.124048 -0.083 0.9345
## ln.vol_cub 0.327657 0.124026 2.642 0.0127 *
## c.con_art5 0.003396 0.053773 0.063 0.9500
## cn_rat -0.066657 0.056740 -1.175 0.2487
## host1:c.mal 1.325556 1.247692 1.062 0.2960
## host1:c.con_art5 0.071246 0.062926 1.132 0.2660
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.635073)
##
## Null deviance: 84.078 on 39 degrees of freedom
## Residual deviance: 52.518 on 32 degrees of freedom
## (2 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
# final model
an3.1 <- glm(rich_art_5 ~ host*c.mal + ln.vol_cub + cn2 + host*c.con_art5,
family=quasipoisson, data=n1)
summary(an3.1)
##
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + ln.vol_cub + cn2 +
## host * c.con_art5, family = quasipoisson, data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3927 -1.0121 -0.2293 0.6199 2.1031
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.17357 1.18643 -1.832 0.07076 .
## host1 -0.19521 0.31487 -0.620 0.53708
## c.mal -0.30162 0.34848 -0.866 0.38940
## ln.vol_cub 0.37460 0.05774 6.487 7.32e-09 ***
## cn2 -0.02960 0.03520 -0.841 0.40301
## c.con_art5 -0.02420 0.02222 -1.089 0.27959
## host1:c.mal 1.10013 0.40981 2.684 0.00887 **
## host1:c.con_art5 0.09438 0.03096 3.049 0.00314 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.346103)
##
## Null deviance: 221.89 on 85 degrees of freedom
## Residual deviance: 110.07 on 78 degrees of freedom
## (60 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an3.f <- glm(rich_art_5 ~ host*c.mal + ln.vol_cub + host*c.con_art5,
family=quasipoisson, data=n1)
summary(an3.f)
##
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + ln.vol_cub + host *
## c.con_art5, family = quasipoisson, data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3865 -1.0153 -0.2121 0.7193 2.0200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.77989 0.63926 -4.349 3.75e-05 ***
## host1 0.02279 0.17607 0.129 0.89730
## c.mal -0.25980 0.33514 -0.775 0.44034
## ln.vol_cub 0.35623 0.05413 6.581 3.51e-09 ***
## c.con_art5 -0.02600 0.02146 -1.212 0.22897
## host1:c.mal 0.98839 0.39786 2.484 0.01492 *
## host1:c.con_art5 0.08287 0.02862 2.896 0.00479 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
##
## Null deviance: 234.49 on 92 degrees of freedom
## Residual deviance: 122.51 on 86 degrees of freedom
## (53 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
# partial r-sq
an3 <- glm(rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1, subset=host=="C")
summary(an3) # R-sq = 53% total
##
## Call:
## glm(formula = rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson,
## data = n1, subset = host == "C")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3858 -1.1404 -0.2420 0.8255 2.0197
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.75955 0.79599 -3.467 0.00098 ***
## c.mal 0.72862 0.22775 3.199 0.00220 **
## ln.vol_cub 0.35643 0.06266 5.688 4.05e-07 ***
## c.con_art5 0.05688 0.01992 2.855 0.00590 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.538553)
##
## Null deviance: 184.732 on 63 degrees of freedom
## Residual deviance: 97.756 on 60 degrees of freedom
## (9 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
n1.p <- n1[complete.cases(n1$rich_art_5, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5), ]
n1.p$rich_art_5 <- sqrt(n1.p$rich_art_5)
an.rich_mal <- lm(rich_art_5 ~ ln.vol_cub + c.con_art5,
data=n1.p, subset=host=="C")
an.rich_vol <- lm(rich_art_5 ~ c.mal + c.con_art5,
data=n1.p, subset=host=="C")
an.rich_con <- lm(rich_art_5 ~ c.mal + ln.vol_cub,
data=n1.p, subset=host=="C")
an.mal <- lm(c.mal ~ ln.vol_cub + c.con_art5,
data=n1.p, subset=host=="C")
an.vol <- lm(ln.vol_cub ~ c.mal + c.con_art5,
data=n1.p, subset=host=="C")
an.con <- lm(c.con_art5 ~ c.mal + ln.vol_cub,
data=n1.p, subset=host=="C")
r.rich_mal <- resid(an.rich_mal)
r.rich_vol <- resid(an.rich_vol)
r.rich_con <- resid(an.rich_con)
r.mal <- resid(an.mal)
r.vol <- resid(an.vol)
r.con <- resid(an.con)
an.p_mal <- lm(r.rich_mal ~ r.mal)
summary(an.p_mal) # r-sq = 9.92 %
##
## Call:
## lm(formula = r.rich_mal ~ r.mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.55867 -0.48803 0.00644 0.56164 1.30531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.816e-17 8.486e-02 0.000 1.00000
## r.mal 7.012e-01 2.488e-01 2.818 0.00648 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6789 on 62 degrees of freedom
## Multiple R-squared: 0.1135, Adjusted R-squared: 0.09924
## F-statistic: 7.941 on 1 and 62 DF, p-value: 0.006476
an.p_vol <- lm(r.rich_vol ~ r.vol)
summary(an.p_vol) # r-sq = 42.3 %
##
## Call:
## lm(formula = r.rich_vol ~ r.vol)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.55867 -0.48803 0.00644 0.56164 1.30531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.469e-17 8.486e-02 0.000 1
## r.vol 3.864e-01 5.626e-02 6.868 3.63e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6789 on 62 degrees of freedom
## Multiple R-squared: 0.4321, Adjusted R-squared: 0.4229
## F-statistic: 47.18 on 1 and 62 DF, p-value: 3.631e-09
an.p_con <- lm(r.rich_con ~ r.con)
summary(an.p_con) # r-sq = 5.84 %
##
## Call:
## lm(formula = r.rich_con ~ r.con)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.55867 -0.48803 0.00644 0.56164 1.30531
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.907e-17 8.486e-02 0.000 1.0000
## r.con 6.205e-02 2.801e-02 2.215 0.0304 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6789 on 62 degrees of freedom
## Multiple R-squared: 0.07335, Adjusted R-squared: 0.0584
## F-statistic: 4.908 on 1 and 62 DF, p-value: 0.03042
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)
an5 <- lm(cn2 ~ c.mal + host + ln.vol_cub , data=n1, subset=no_tim > 0)
summary(an5)
##
## Call:
## lm(formula = cn2 ~ c.mal + host + ln.vol_cub, data = n1, subset = no_tim >
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3563 -1.3177 0.0028 1.1640 4.5875
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.94660 1.90134 10.491 <2e-16 ***
## c.mal -1.31193 0.54549 -2.405 0.0184 *
## host1 7.56694 0.46137 16.401 <2e-16 ***
## ln.vol_cub 0.05646 0.15487 0.365 0.7164
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.81 on 82 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.8068, Adjusted R-squared: 0.7997
## F-statistic: 114.1 on 3 and 82 DF, p-value: < 2.2e-16
an6 <- lm(cn2 ~ pC, data=n1)
an7 <- lm(cn2 ~ pN, data=n1)
summary(an6)
##
## Call:
## lm(formula = cn2 ~ pC, data = n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.898 -3.701 -1.085 3.810 14.042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.0853 3.2998 8.208 1.59e-13 ***
## pC -0.2731 0.3185 -0.858 0.393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.656 on 135 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.005417, Adjusted R-squared: -0.00195
## F-statistic: 0.7353 on 1 and 135 DF, p-value: 0.3927
summary(an7)
##
## Call:
## lm(formula = cn2 ~ pN, data = n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0101 -1.7586 -0.3754 1.3619 18.4994
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.854 1.030 38.70 <2e-16 ***
## pN -35.430 2.279 -15.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.795 on 135 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.6417, Adjusted R-squared: 0.639
## F-statistic: 241.7 on 1 and 135 DF, p-value: < 2.2e-16
contrasts(n1$host) <- c(0, 1)
an5 <- lm(pC ~ c.mal + host, data=n1, subset=no_tim>0)
summary(an5)
##
## Call:
## lm(formula = pC ~ c.mal + host, data = n1, subset = no_tim >
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0288 -0.5477 -0.0229 0.5032 11.0203
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.1707 0.2980 34.132 <2e-16 ***
## c.mal 0.2730 0.4463 0.612 0.542
## host1 0.2438 0.3606 0.676 0.501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.482 on 83 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.01341, Adjusted R-squared: -0.01036
## F-statistic: 0.564 on 2 and 83 DF, p-value: 0.5711
contrasts(n1$host) <- c(0, 1)
an5 <- lm(pN ~ c.mal + host + no_herb , data=n1, subset=no_tim>0)
summary(an5)
##
## Call:
## lm(formula = pN ~ c.mal + host + no_herb, data = n1, subset = no_tim >
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15346 -0.03789 0.00174 0.02453 0.58215
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3535096 0.0202036 17.497 < 2e-16 ***
## c.mal 0.0380642 0.0251261 1.515 0.134
## host1 0.1418512 0.0199923 7.095 4.19e-10 ***
## no_herb 0.0008782 0.0008526 1.030 0.306
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08202 on 82 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.4475, Adjusted R-squared: 0.4273
## F-statistic: 22.14 on 3 and 82 DF, p-value: 1.35e-10
# final model
an5.1 <- lm(cn2 ~ c.mal + host + ln.vol_cub + no_art_5, data=n1, subset=no_tim > 0)
summary(an5.1)
##
## Call:
## lm(formula = cn2 ~ c.mal + host + ln.vol_cub + no_art_5, data = n1,
## subset = no_tim > 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3052 -1.1541 -0.0549 1.0570 4.3339
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.37517 1.91214 13.794 <2e-16 ***
## c.mal -1.13413 0.55220 -2.054 0.0432 *
## host1 -7.71009 0.46613 -16.541 <2e-16 ***
## ln.vol_cub 0.19395 0.17646 1.099 0.2750
## no_art_5 -0.04080 0.02585 -1.579 0.1183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.794 on 81 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.8125, Adjusted R-squared: 0.8033
## F-statistic: 87.76 on 4 and 81 DF, p-value: < 2.2e-16
an5.2 <- lm(cn2 ~ c.mal + host + no_art_5, data=n1, subset=no_tim > 0)
summary(an5.2)# R-sq = 80.2 %
##
## Call:
## lm(formula = cn2 ~ c.mal + host + no_art_5, data = n1, subset = no_tim >
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.377 -1.086 -0.082 1.093 4.516
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.42373 0.42761 66.471 <2e-16 ***
## c.mal -1.17752 0.55149 -2.135 0.0357 *
## host1 -7.53107 0.43730 -17.222 <2e-16 ***
## no_art_5 -0.02678 0.02251 -1.190 0.2375
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.796 on 82 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.8097, Adjusted R-squared: 0.8028
## F-statistic: 116.3 on 3 and 82 DF, p-value: < 2.2e-16
an5.f <- lm(cn2 ~ c.mal + host, data=n1, subset=no_tim > 0)
summary(an5.f)
##
## Call:
## lm(formula = cn2 ~ c.mal + host, data = n1, subset = no_tim >
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.377 -1.280 -0.009 1.170 4.624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.1514 0.3621 77.739 <2e-16 ***
## c.mal -1.3051 0.5423 -2.407 0.0183 *
## host1 -7.5170 0.4382 -17.153 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.8 on 83 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.8064, Adjusted R-squared: 0.8018
## F-statistic: 172.9 on 2 and 83 DF, p-value: < 2.2e-16
## neg freq dep mal
an5.f <- lm(cn2 ~ c.nf_mal_k0.5 + host, data=n1, subset=no_tim > 0)
summary(an5.f)
##
## Call:
## lm(formula = cn2 ~ c.nf_mal_k0.5 + host, data = n1, subset = no_tim >
## 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1498 -1.3665 -0.0076 1.2264 4.2070
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.2882 0.3566 79.330 <2e-16 ***
## c.nf_mal_k0.5 -1.1872 0.5567 -2.133 0.0359 *
## host1 -7.7235 0.4270 -18.086 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.813 on 83 degrees of freedom
## (8 observations deleted due to missingness)
## Multiple R-squared: 0.8037, Adjusted R-squared: 0.799
## F-statistic: 169.9 on 2 and 83 DF, p-value: < 2.2e-16
##
# partial r-sq
p.n1 <- n1[complete.cases(n1$host, n1$c.mal, n1$cn2), ]
p.n1$c.host <- ifelse(p.n1$host == "A", -.5, .5)
an.mal <- lm(c.mal ~ c.host, data=p.n1)
an.host <- lm(c.host ~ c.mal, data=p.n1)
an.cn_h <- lm(cn2 ~ c.host, data=p.n1)
an.cn_m <- lm(cn2 ~ c.mal, data=p.n1)
r.mal <- resid(an.mal)
r.host <- resid(an.host)
r.cn_h <- resid(an.cn_h)
r.cn_m <- resid(an.cn_m)
an.p_mal <- lm(r.cn_h ~ r.mal)
summary(an.p_mal) # r-sq = 5.41%
##
## Call:
## lm(formula = r.cn_h ~ r.mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.377 -1.280 -0.009 1.170 4.624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.202e-17 1.930e-01 0.000 1.0000
## r.mal -1.305e+00 5.391e-01 -2.421 0.0176 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.79 on 84 degrees of freedom
## Multiple R-squared: 0.06523, Adjusted R-squared: 0.0541
## F-statistic: 5.861 on 1 and 84 DF, p-value: 0.01763
an.p_host <- lm(r.cn_m ~ r.host)
summary(an.p_host) # r-sq = 77.7%
##
## Call:
## lm(formula = r.cn_m ~ r.host)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.377 -1.280 -0.009 1.170 4.624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.551e-16 1.930e-01 0.00 1
## r.host -7.517e+00 4.356e-01 -17.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.79 on 84 degrees of freedom
## Multiple R-squared: 0.78, Adjusted R-squared: 0.7773
## F-statistic: 297.8 on 1 and 84 DF, p-value: < 2.2e-16
# on biomass
# changed starts from 0
an2 <- glm(biomass ~ mal, family=Gamma, start=c(0.0001,0.0001),
data=n1[n1$biomass > 0,])
summary(an2)
##
## Call:
## glm(formula = biomass ~ mal, family = Gamma, data = n1[n1$biomass >
## 0, ], start = c(1e-04, 1e-04))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4097 -1.0986 -0.4272 0.2647 2.0929
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.295 1.240 5.885 7.47e-08 ***
## mal -3.431 1.733 -1.980 0.0509 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 1.196149)
##
## Null deviance: 111.56 on 87 degrees of freedom
## Residual deviance: 106.84 on 86 degrees of freedom
## (41 observations deleted due to missingness)
## AIC: -120.38
##
## Number of Fisher Scoring iterations: 21
## connectivity on maladaptation?
mal_nomal <- cbind(n1_occ$no_mal, n1_occ$no_adp)
grn_str <- cbind(n1_occ$no_grn_tim, n1_occ$no_str_tim)
#an1 <- glm(tims ~ n1_occ$con_tim, family=binomial) # no effect of connectivity, !!what is tims
an2 <- glm(mal ~ con_tim, data=n1_occ)
summary(an2)
##
## Call:
## glm(formula = mal ~ con_tim, data = n1_occ)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.49278 -0.43417 0.02203 0.29194 0.56315
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.49289 0.05898 8.357 6.96e-13 ***
## con_tim -0.01668 0.02085 -0.800 0.426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1430296)
##
## Null deviance: 13.107 on 92 degrees of freedom
## Residual deviance: 13.016 on 91 degrees of freedom
## AIC: 87.043
##
## Number of Fisher Scoring iterations: 2
an3 <- glm(mal_nomal ~ no_art_5, family=binomial, data=n1_occ, subset=no_tim>0)
summary(an3)
##
## Call:
## glm(formula = mal_nomal ~ no_art_5, family = binomial, data = n1_occ,
## subset = no_tim > 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5271 -1.0306 -0.3327 0.8986 2.5854
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.338683 0.166728 -2.031 0.0422 *
## no_art_5 -0.002842 0.010716 -0.265 0.7909
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 133.15 on 92 degrees of freedom
## Residual deviance: 133.08 on 91 degrees of freedom
## AIC: 239.92
##
## Number of Fisher Scoring iterations: 3
an4 <- lm(mal ~ no_art_5 + ln.vol_cub + host + c.con_art5, data=n1)
summary(an4)
##
## Call:
## lm(formula = mal ~ no_art_5 + ln.vol_cub + host + c.con_art5,
## data = n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.55503 -0.27132 -0.03773 0.27866 0.76125
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.515242 0.330751 1.558 0.12287
## no_art_5 0.009593 0.004962 1.933 0.05643 .
## ln.vol_cub -0.027387 0.030703 -0.892 0.37482
## host1 0.248194 0.094221 2.634 0.00996 **
## c.con_art5 -0.008293 0.009454 -0.877 0.38281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3568 on 88 degrees of freedom
## (53 observations deleted due to missingness)
## Multiple R-squared: 0.1451, Adjusted R-squared: 0.1062
## F-statistic: 3.734 on 4 and 88 DF, p-value: 0.007463
# plant growth
hist(n1$avg_growth)
an1 <- lm(avg_growth ~ con_art5 + host + no_art_5, data=n1)
summary(an1)
##
## Call:
## lm(formula = avg_growth ~ con_art5 + host + no_art_5, data = n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.527 -12.393 -1.892 10.114 45.407
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 81.8856 4.0065 20.438 <2e-16 ***
## con_art5 -0.9006 0.3759 -2.396 0.0183 *
## host1 -40.4208 3.9231 -10.303 <2e-16 ***
## no_art_5 -0.5023 0.2250 -2.232 0.0277 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.14 on 105 degrees of freedom
## (37 observations deleted due to missingness)
## Multiple R-squared: 0.5578, Adjusted R-squared: 0.5452
## F-statistic: 44.15 on 3 and 105 DF, p-value: < 2.2e-16
hist(n1$avg_no_flwrs)
n1$no_not_herb_5 <- n1$no_art_5 - n1$no_herb_5
an2 <- glm(avg_no_flwrs ~ no_herb_5 + ln.vol_cub, data=n1, family=quasipoisson)
summary(an2)
##
## Call:
## glm(formula = avg_no_flwrs ~ no_herb_5 + ln.vol_cub, family = quasipoisson,
## data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8752 -1.0093 -0.5487 0.2925 2.5442
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.53768 1.22170 -2.896 0.00508 **
## no_herb_5 0.01450 0.01187 1.221 0.22618
## ln.vol_cub 0.29910 0.10407 2.874 0.00540 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.286064)
##
## Null deviance: 109.804 on 70 degrees of freedom
## Residual deviance: 83.553 on 68 degrees of freedom
## (75 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 6
contrasts(n1$host) <- c(-0.5, 0.5)
contrasts(n1$host) <- c(0, 1)
contrasts(n1$host) <- c(1, 0)
an.s1 <- lm(mas ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim + c.con_art5, data=n1[(n1$no_art_5 > 2 & n1$no_tim > 0 & n1$c.mal < 0.6), ])
summary(an.s1)
##
## Call:
## lm(formula = mas ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim + c.con_art5, data = n1[(n1$no_art_5 > 2 & n1$no_tim >
## 0 & n1$c.mal < 0.6), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60977 -0.18461 -0.01867 0.15665 0.64327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.681677 0.053586 -12.721 < 2e-16 ***
## host1 0.214007 0.103443 2.069 0.04242 *
## c.mal 0.150830 0.129742 1.163 0.24914
## c.no_art_5 -0.030005 0.004684 -6.406 1.72e-08 ***
## c.no_tim -0.014150 0.009853 -1.436 0.15560
## c.con_art5 0.027374 0.008206 3.336 0.00139 **
## host1:c.mal 0.116187 0.228700 0.508 0.61310
## host1:c.no_art_5 -0.021955 0.010567 -2.078 0.04157 *
## c.mal:c.no_tim 0.101103 0.036701 2.755 0.00755 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2837 on 67 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.5893, Adjusted R-squared: 0.5403
## F-statistic: 12.02 on 8 and 67 DF, p-value: 1.833e-10
an.sf <- lm(mas ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim + c.con_art5, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an.sf)
##
## Call:
## lm(formula = mas ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim + c.con_art5, data = n1[n1$no_art_5 > 2 & n1$no_tim >
## 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60977 -0.18461 -0.01867 0.15665 0.64327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.681677 0.053586 -12.721 < 2e-16 ***
## host1 0.214007 0.103443 2.069 0.04242 *
## c.mal 0.150830 0.129742 1.163 0.24914
## c.no_art_5 -0.030005 0.004684 -6.406 1.72e-08 ***
## c.no_tim -0.014150 0.009853 -1.436 0.15560
## c.con_art5 0.027374 0.008206 3.336 0.00139 **
## host1:c.mal 0.116187 0.228700 0.508 0.61310
## host1:c.no_art_5 -0.021955 0.010567 -2.078 0.04157 *
## c.mal:c.no_tim 0.101103 0.036701 2.755 0.00755 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2837 on 67 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.5893, Adjusted R-squared: 0.5403
## F-statistic: 12.02 on 8 and 67 DF, p-value: 1.833e-10
## neg freq dep mal
#an.sf <- lm(mas ~ host*c.nf_mal_k2.0 + host*c.no_art_5 + c.nf_mal_k2.0*c.no_tim + c.con_art5, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
#summary(an.sf)
###
an1 <- lm(mai1 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an1)
##
## Call:
## lm(formula = mai1 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.50939 -0.40094 -0.02639 0.31742 2.32790
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.488093 0.115158 4.238 6.94e-05 ***
## host1 -0.507877 0.212174 -2.394 0.01945 *
## c.mal 0.054880 0.283720 0.193 0.84720
## c.no_art_5 0.052519 0.010370 5.065 3.34e-06 ***
## c.no_tim 0.004581 0.021266 0.215 0.83010
## host1:c.mal -1.054458 0.502242 -2.100 0.03949 *
## host1:c.no_art_5 0.044190 0.023411 1.888 0.06335 .
## c.mal:c.no_tim -0.276703 0.080579 -3.434 0.00102 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6286 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4525, Adjusted R-squared: 0.3961
## F-statistic: 8.028 on 7 and 68 DF, p-value: 4.402e-07
an20 <- lm(mai20 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an20)
##
## Call:
## lm(formula = mai20 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3788 -0.4707 -0.0875 0.3971 3.4183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.06086 0.17739 -0.343 0.732602
## host1 1.22800 0.32684 3.757 0.000359 ***
## c.mal 0.03294 0.43705 0.075 0.940145
## c.no_art_5 -0.01297 0.01597 -0.812 0.419549
## c.no_tim 0.02004 0.03276 0.612 0.542669
## host1:c.mal 2.45150 0.77366 3.169 0.002295 **
## host1:c.no_art_5 -0.06586 0.03606 -1.826 0.072183 .
## c.mal:c.no_tim 0.33483 0.12413 2.698 0.008801 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9683 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2636, Adjusted R-squared: 0.1878
## F-statistic: 3.477 on 7 and 68 DF, p-value: 0.003047
an3 <- lm(mai3 ~ host*c.mal+ host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an3)
##
## Call:
## lm(formula = mai3 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.99069 -0.30888 -0.04333 0.22945 1.86612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.430309 0.088111 4.884 6.62e-06 ***
## host1 -0.325153 0.162341 -2.003 0.049177 *
## c.mal 0.052570 0.217083 0.242 0.809381
## c.no_art_5 0.045625 0.007934 5.750 2.31e-07 ***
## c.no_tim 0.006208 0.016271 0.382 0.703982
## host1:c.mal -0.685410 0.384281 -1.784 0.078948 .
## host1:c.no_art_5 0.032606 0.017913 1.820 0.073120 .
## c.mal:c.no_tim -0.212331 0.061654 -3.444 0.000986 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4809 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4951, Adjusted R-squared: 0.4432
## F-statistic: 9.528 on 7 and 68 DF, p-value: 3.449e-08
an5 <- lm(mai5 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an5)
##
## Call:
## lm(formula = mai5 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5089 -0.2210 -0.0606 0.1439 1.4043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.372524 0.063327 5.883 1.36e-07 ***
## host1 -0.142429 0.116677 -1.221 0.22641
## c.mal 0.050260 0.156021 0.322 0.74834
## c.no_art_5 0.038731 0.005703 6.792 3.36e-09 ***
## c.no_tim 0.007836 0.011694 0.670 0.50508
## host1:c.mal -0.316361 0.276189 -1.145 0.25604
## host1:c.no_art_5 0.021021 0.012874 1.633 0.10713
## c.mal:c.no_tim -0.147959 0.044311 -3.339 0.00137 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3457 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.558, Adjusted R-squared: 0.5125
## F-statistic: 12.27 on 7 and 68 DF, p-value: 4.965e-10
an7 <- lm(mai7 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an7)
##
## Call:
## lm(formula = mai7 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38411 -0.15051 -0.05102 0.08080 0.94254
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.314740 0.044738 7.035 1.23e-09 ***
## host1 0.040294 0.082429 0.489 0.62653
## c.mal 0.047951 0.110224 0.435 0.66492
## c.no_art_5 0.031837 0.004029 7.903 3.29e-11 ***
## c.no_tim 0.009464 0.008262 1.145 0.25602
## host1:c.mal 0.052687 0.195119 0.270 0.78796
## host1:c.no_art_5 0.009436 0.009095 1.037 0.30318
## c.mal:c.no_tim -0.083587 0.031305 -2.670 0.00948 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2442 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.623, Adjusted R-squared: 0.5842
## F-statistic: 16.05 on 7 and 68 DF, p-value: 2.902e-12
an8 <- lm(mai8 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an8)
##
## Call:
## lm(formula = mai8 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33478 -0.16438 -0.03142 0.10841 0.74215
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.285848 0.040637 7.034 1.23e-09 ***
## host1 0.131656 0.074872 1.758 0.0832 .
## c.mal 0.046796 0.100119 0.467 0.6417
## c.no_art_5 0.028390 0.003659 7.758 6.02e-11 ***
## c.no_tim 0.010277 0.007504 1.370 0.1753
## host1:c.mal 0.237211 0.177231 1.338 0.1852
## host1:c.no_art_5 0.003644 0.008261 0.441 0.6606
## c.mal:c.no_tim -0.051401 0.028435 -1.808 0.0751 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2218 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6229, Adjusted R-squared: 0.584
## F-statistic: 16.04 on 7 and 68 DF, p-value: 2.933e-12
an9 <- lm(mai9 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an9)
##
## Call:
## lm(formula = mai9 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38322 -0.15543 -0.04908 0.11741 0.70110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.256956 0.041681 6.165 4.37e-08 ***
## host1 0.223018 0.076796 2.904 0.00496 **
## c.mal 0.045641 0.102692 0.444 0.65813
## c.no_art_5 0.024943 0.003753 6.646 6.13e-09 ***
## c.no_tim 0.011091 0.007697 1.441 0.15418
## host1:c.mal 0.421735 0.181785 2.320 0.02335 *
## host1:c.no_art_5 -0.002149 0.008474 -0.254 0.80060
## c.mal:c.no_tim -0.019215 0.029165 -0.659 0.51222
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2275 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.5745, Adjusted R-squared: 0.5307
## F-statistic: 13.12 on 7 and 68 DF, p-value: 1.461e-10
an10 <- lm(mai10 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an10)
##
## Call:
## lm(formula = mai10 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45265 -0.16372 -0.05317 0.14843 0.82478
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.228063 0.047533 4.798 9.13e-06 ***
## host1 0.314380 0.087578 3.590 0.00062 ***
## c.mal 0.044486 0.117110 0.380 0.70523
## c.no_art_5 0.021497 0.004280 5.022 3.93e-06 ***
## c.no_tim 0.011905 0.008778 1.356 0.17950
## host1:c.mal 0.606259 0.207308 2.924 0.00468 **
## host1:c.no_art_5 -0.007941 0.009663 -0.822 0.41409
## c.mal:c.no_tim 0.012971 0.033260 0.390 0.69778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2595 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4891, Adjusted R-squared: 0.4365
## F-statistic: 9.3 on 7 and 68 DF, p-value: 5.018e-08
an18 <- lm(mai18 ~ host*c.mal + host*c.no_art_5 + c.mal*c.no_tim,
data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
summary(an18)
##
## Call:
## lm(formula = mai18 ~ host * c.mal + host * c.no_art_5 + c.mal *
## c.no_tim, data = n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.98226 -0.39681 -0.08116 0.33437 2.89961
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.003074 0.148785 -0.021 0.983578
## host1 1.045274 0.274130 3.813 0.000298 ***
## c.mal 0.035248 0.366569 0.096 0.923679
## c.no_art_5 -0.006079 0.013398 -0.454 0.651466
## c.no_tim 0.018416 0.027476 0.670 0.504960
## host1:c.mal 2.082452 0.648900 3.209 0.002032 **
## host1:c.no_art_5 -0.054280 0.030247 -1.795 0.077173 .
## c.mal:c.no_tim 0.270458 0.104109 2.598 0.011491 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8121 on 68 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2619, Adjusted R-squared: 0.1859
## F-statistic: 3.446 on 7 and 68 DF, p-value: 0.00325
# predictions on Adenostoma
pred.slope <- predict(an.sf, newdata=data.frame(host="A", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE) # changed an.s for an.sf
pred.int <- predict(an1, newdata=data.frame(host="A", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE)
pred.int20 <- predict(an20, newdata=data.frame(host="A", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE)
# predictions on Ceanothus - do not plot as not significant
# pred.slope <- predict(an.sf, newdata=data.frame(host="C", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE) # changed an.s for an.s1
# pred.int <- predict(an1, newdata=data.frame(host="C", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0, c.con_art5=0), se.fit=TRUE)
# pred.int20 <- predict(an20, newdata=data.frame(host="C", c.mal=c(-.5, .5), c.no_art_5=0, c.no_tim=0), c.con_art5=0, se.fit=TRUE)
# partial r-sq
an.s2 <- lm(mas ~ c.mal + c.no_art_5, data=n1[(n1$no_art_5 > 2 & n1$no_tim > 0 & n1$host == "A"), ])
summary(an.s2) # R-sq = 18.4 %
##
## Call:
## lm(formula = mas ~ c.mal + c.no_art_5, data = n1[(n1$no_art_5 >
## 2 & n1$no_tim > 0 & n1$host == "A"), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62088 -0.20001 0.03397 0.15980 0.85752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.40958 0.10117 -4.048 0.000536 ***
## c.mal 0.18082 0.21821 0.829 0.416207
## c.no_art_5 -0.04825 0.01169 -4.127 0.000443 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3531 on 22 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.4409, Adjusted R-squared: 0.3901
## F-statistic: 8.675 on 2 and 22 DF, p-value: 0.001668
n1.p <- n1[complete.cases(n1$mas, n1$c.mal,n1$c.no_art_5), ]
# to drop levels properly, keep vars of interest
n1.p <- n1.p[,c("c.no_art_5", "c.mal", "no_tim", "host", "no_art_5", "mas")]
an.mas.art <- lm(mas ~ c.no_art_5, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
an.mas.mal <- lm(mas ~ c.mal, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
an.mal.art <- lm(c.mal ~ c.no_art_5, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
an.art.mal <- lm(c.no_art_5 ~ c.mal, data=n1.p[(n1.p$no_art_5 > 2 & n1.p$no_tim > 0 & n1.p$host == "A"), ])
r.mas <- resid(an.mas.art)
r.mal <- resid(an.mal.art)
r.art <- resid(an.art.mal)
an.p_mal <- lm(r.mas ~ r.mal)
summary(an.p_mal) # r-sq = 0.4% on Adenostoma
##
## Call:
## lm(formula = r.mas ~ r.mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62088 -0.20001 0.03397 0.15980 0.85752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.353e-18 6.906e-02 0.000 1.000
## r.mal 1.808e-01 2.134e-01 0.847 0.406
##
## Residual standard error: 0.3453 on 23 degrees of freedom
## Multiple R-squared: 0.03027, Adjusted R-squared: -0.0119
## F-statistic: 0.7179 on 1 and 23 DF, p-value: 0.4056
an.p_art <- lm(r.mal~ r.art)
summary(an.p_art) # r-sq = 7% on Adenostoma
##
## Call:
## lm(formula = r.mal ~ r.art)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2155 -0.2155 -0.2155 0.1664 0.6756
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.051e-17 6.370e-02 0.000 1.000
## r.art -1.768e-02 1.055e-02 -1.676 0.107
##
## Residual standard error: 0.3185 on 23 degrees of freedom
## Multiple R-squared: 0.1089, Adjusted R-squared: 0.07012
## F-statistic: 2.81 on 1 and 23 DF, p-value: 0.1072
pdf(file="figures/mas_curves_by_mal4A.pdf")
plot(x=c(0, 19), y=c(-1, 3), type="n", las=1,
xlab="rank log biomass", ylab="log abundance", axes=FALSE)
axis(1, at=c(0, 4, 9, 14, 19), labels=c(1, 5, 10, 15, 20))
axis(2, at=seq(-1, 3, by=.5),
labels=c("-1", "", "0", "", "1", "", "2", "", "3"), las=1)
# lines(x=c(0, 19), y=c(pred.int$fit[1]-pred.int$se.fit[1],
# pred.int20$fit[1]+pred.int$se.fit[1]),
# col="blue", lty=2)
# lines(x=c(0, 19), y=c(pred.int$fit[1]+pred.int$se.fit[1],
# pred.int20$fit[1]-pred.int$se.fit[1]),
# col="blue", lty=2)
polygon(x=c(0, 0, 19, 19), y=c(pred.int$fit[1] - pred.int$se.fit[1],
pred.int$fit[1] + pred.int$se.fit[1],
pred.int20$fit[1] - pred.int20$se.fit[1],
pred.int20$fit[1] + pred.int20$se.fit[1]),
lwd=.5, col="blue", density=150)
lines(x=c(0, 19), y=c(pred.int$fit[1], pred.int20$fit[1]),
col="blue", lwd=2)
# lines(x=c(0, 19), y=c(pred.int$fit[2]-pred.int$se.fit[2],
# pred.int20$fit[2]+pred.int$se.fit[2]),
# col="red", lty=2)
# lines(x=c(0, 19), y=c(pred.int$fit[2]+pred.int$se.fit[2],
# pred.int20$fit[2]-pred.int$se.fit[2]),
# col="red", lty=2)
polygon(x=c(0, 0, 19, 19), y=c(pred.int$fit[2] - pred.int$se.fit[2],
pred.int$fit[2] + pred.int$se.fit[2],
pred.int20$fit[2] - pred.int20$se.fit[2],
pred.int20$fit[2] + pred.int20$se.fit[2]),
lwd=.5, col="red", density=150)
lines(x=c(0, 19), y=c(pred.int$fit[2], pred.int20$fit[2]),
col="red", lwd=2)
abline(v=c(2, 8), lty=2, lwd=1, col="darkgreen")
legend(3, 2.5, legend=c("maladapted", "well adapted"), fill=c("red", "blue"),
bty="n", cex=.9)
arrows(x0=c(2, 8), x1=c(0, 11), y0=c(0.075, 0.5), y1=c(0.075, 0.5), length=.1)
text(x=c(1, 9.25), y=c(0.2, .6), labels=c("p < 0.05", "p < 0.05"), cex=0.85)
dev.off()
## quartz_off_screen
## 2
contrasts(n1$host) <- c(0, 1)
an1 <- lm(mas ~ no_art_5 + c.mal*host, data=n1[n1$mas<0.6 & n1$no_art_5 > 5, ])
summary(an1)
##
## Call:
## lm(formula = mas ~ no_art_5 + c.mal * host, data = n1[n1$mas <
## 0.6 & n1$no_art_5 > 5, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.54968 -0.17396 -0.06562 0.16417 1.00561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.29509 0.10174 -2.900 0.00538 **
## no_art_5 -0.02844 0.00529 -5.376 1.67e-06 ***
## c.mal 0.07591 0.19195 0.395 0.69406
## host1 -0.28841 0.09800 -2.943 0.00478 **
## c.mal:host1 0.20087 0.26738 0.751 0.45577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3186 on 54 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.4593, Adjusted R-squared: 0.4193
## F-statistic: 11.47 on 4 and 54 DF, p-value: 8.242e-07
an1 <- lm(mas ~ no_art_5, data=n1)
summary(an1)
##
## Call:
## lm(formula = mas ~ no_art_5, data = n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63760 -0.20930 -0.06227 0.17180 1.12303
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.37902 0.04887 -7.755 6.00e-12 ***
## no_art_5 -0.03222 0.00368 -8.755 3.75e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3261 on 105 degrees of freedom
## (39 observations deleted due to missingness)
## Multiple R-squared: 0.422, Adjusted R-squared: 0.4165
## F-statistic: 76.66 on 1 and 105 DF, p-value: 3.748e-14
plot(n1$mas[n1$no_art_5 > 5] ~ n1$no_art_5[n1$no_art_5 > 5])
nfds.stats <- as.data.frame(matrix(nrow = 7, ncol=12))
# Timema abundance
for(i in 1:6) {
an.tim <- glm(n1$no_tim ~ n1[[189 + i]] + n1$ln.vol_cub + n1$con_tim,
family=quasipoisson, data=n1,
subset=no_tim > 4)
nfds.stats[1, c(2 * i - 1, 2 * i)] <- summary(an.tim)$coefficients[2, c(1, 4)]
}
# arthropod abundance
for(i in 1:6) {
an.art <- glm(no_art_5 ~ n1[[189+i]] + n1$ln.vol_cub + n1$c.con_art5*n1$host,
family=quasipoisson, data=n1, subset = no_tim > 0)
nfds.stats[2, c(2 * i - 1, 2 * i)] <- summary(an.art)$coefficients[2, c(1, 4)]
}
# arthropod richness
contrasts(n1$host) <- c(1, 0)
for(i in 1:6) {
an.rich <- glm(n1$rich_art_5 ~ n1$host * n1[[189+i]] + n1$ln.vol_cub + n1$cn2 + n1$host * n1$c.con_art5,
family=quasipoisson, data=n1)
nfds.stats[3:4, c(2 * i - 1, 2 * i)] <- summary(an.rich)$coefficients[c(3, 7), c(1, 4)]
}
# Mass-abundance
contrasts(n1$host) <- c(0, 1)
for(i in 1:6) {
an.mas <- lm(n1$mas ~ n1$host * n1[[189 + i]] + n1$host * n1$c.no_art_5 + n1[[189 + i]] * n1$c.no_tim + n1$c.con_art5, data=n1[n1$no_art_5 > 2 & n1$no_tim > 0, ])
nfds.stats[5:6, c(2 * i - 1, 2 * i)] <- summary(an.mas)$coefficients[c(3, 7), c(1, 4)]
}
# carbon - nitrogen
for(i in 1:6) {
an.cn <- lm(n1$cn2 ~ n1[[189 + i]] + n1$host + n1$ln.vol_cub + n1$no_art_5, data=n1, subset=no_tim > 0)
nfds.stats[7, c(2 * i - 1, 2 * i)] <- summary(an.cn)$coefficients[2, c(1, 4)]
}
nfds.stats$term <- c("mal", "mal", "mal", "mal * host", "mal", "mal * host", "mal")
colnames(nfds.stats) <- rep(c("b", "p"), 6)
write.csv(nfds.stats, file = "data/nfds_stats.csv")
Figures
jpeg(file="figures/n1_tim_mal.jpg",
width=750, height=750)
par(mar=c(5, 6, 1, 1))
plot(y=jitter(n1$no_tim[n1$no_tim > 4], 2),
x=jitter(n1$mal[n1$no_tim > 4], 10),
xlab="maladaptation", ylab="",
ylim=c(3, 17),
las=1, pch=16,
cex=1.5, cex.lab=1.5, cex.axis=1.5)
title(ylab=expression(paste(italic("T. cristinae"), " abundance")), cex.lab=1.5, line=3.5)
axis.break(axis=2, breakpos=3.5)
dev.off()
## quartz_off_screen
## 2
an2.f <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5*host,
family=quasipoisson, data=n1, subset = no_tim > 0)
summary(an2.f)
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5 * host,
## family = quasipoisson, data = n1, subset = no_tim > 0)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6751 -1.6447 -0.7855 0.8684 7.8261
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.594e+00 9.414e-01 -2.755 0.00714 **
## c.mal 4.999e-01 2.369e-01 2.110 0.03770 *
## ln.vol_cub 4.316e-01 8.042e-02 5.367 6.55e-07 ***
## c.con_art5 -8.682e-05 2.590e-02 -0.003 0.99733
## host1 -4.782e-01 2.173e-01 -2.200 0.03042 *
## c.con_art5:host1 6.740e-02 3.701e-02 1.821 0.07201 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 5.346001)
##
## Null deviance: 652.94 on 92 degrees of freedom
## Residual deviance: 407.37 on 87 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
an2.f1 <- glm(no_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1, subset = no_tim > 0 & host== "C")
summary(an2.f1)
##
## Call:
## glm(formula = no_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson,
## data = n1, subset = no_tim > 0 & host == "C")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.6270 -1.5586 -0.7073 0.7373 7.9524
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.05242 1.26538 -2.412 0.0189 *
## c.mal 0.71362 0.34010 2.098 0.0401 *
## ln.vol_cub 0.42820 0.09911 4.320 5.95e-05 ***
## c.con_art5 0.07210 0.02852 2.528 0.0141 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 6.078791)
##
## Null deviance: 512.65 on 63 degrees of freedom
## Residual deviance: 311.59 on 60 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
n1.2 <- n1[complete.cases(n1$no_art_5, n1$host, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5, n1$no_tim), ]
res.art <- resid(glm(no_art_5 ~ ln.vol_cub + c.con_art5*host,
family=quasipoisson, data=n1.2, subset = no_tim > 0 ), 'deviance')
res.mal <- resid(lm(c.mal ~ ln.vol_cub + c.con_art5*host,
data=n1.2, subset = no_tim > 0), "pearson")
an.res <- lm(res.art ~ res.mal); summary(an.res)
##
## Call:
## lm(formula = res.art ~ res.mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4703 -1.3778 -0.3566 1.1054 7.9407
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2809 0.2192 -1.281 0.2033
## res.mal 1.2042 0.6201 1.942 0.0552 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.114 on 91 degrees of freedom
## Multiple R-squared: 0.0398, Adjusted R-squared: 0.02925
## F-statistic: 3.772 on 1 and 91 DF, p-value: 0.05522
nd=data.frame(res.mal=seq(min(res.mal), max(res.mal), len=50))
preds <- predict(an.res, newdata=nd, type = "response", se.fit=TRUE)
nd$fit <- preds$fit
nd$lower <- preds$fit - preds$se.fit
nd$upper <- preds$fit + preds$se.fit
jpeg(file="figures/n1_art_mal.jpg",
width=750, height=750)
par(mar=c(6, 6, 1, 1))
plot(y=res.art, x=res.mal,
xlab="", ylab="",
las=1, pch=16,
cex=2, cex.lab=2, cex.axis=2)
title(ylab="arthropod abundance (residuals)", cex.lab=3, line=3.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4.5)
#axis.break(axis=2, breakpos=3.5)
polygon(x = c(nd$res.mal, rev(nd$res.mal)),
y = c(nd$upper, rev(nd$lower)),
col="grey", density=100, lwd=.5)
lines(nd$res.mal, nd$fit, lwd=2)
dev.off()
## quartz_off_screen
## 2
#lines(nd$res.mal, nd$upper, lty=2)
#lines(nd$res.mal, nd$lower, lty=2)
an3.f <- glm(rich_art_5 ~ host*c.mal + ln.vol_cub + host*c.con_art5,
family=quasipoisson, data=n1)
summary(an3.f)
##
## Call:
## glm(formula = rich_art_5 ~ host * c.mal + ln.vol_cub + host *
## c.con_art5, family = quasipoisson, data = n1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3865 -1.0153 -0.2121 0.7193 2.0200
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.77989 0.63926 -4.349 3.75e-05 ***
## host1 0.02279 0.17607 0.129 0.89730
## c.mal -0.25980 0.33514 -0.775 0.44034
## ln.vol_cub 0.35623 0.05413 6.581 3.51e-09 ***
## c.con_art5 -0.02600 0.02146 -1.212 0.22897
## host1:c.mal 0.98839 0.39786 2.484 0.01492 *
## host1:c.con_art5 0.08287 0.02862 2.896 0.00479 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.359015)
##
## Null deviance: 234.49 on 92 degrees of freedom
## Residual deviance: 122.51 on 86 degrees of freedom
## (53 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
n1.2 <- n1[complete.cases(n1$rich_art_5, n1$c.mal, n1$ln.vol_cub, n1$c.con_art5), ]
# Adenostoma
an3.fA <- glm(rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1, subset=host=="A")
summary(an3.fA)
##
## Call:
## glm(formula = rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson,
## data = n1, subset = host == "A")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6979 -0.7854 -0.1862 0.5058 1.7825
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.76774 1.34874 -2.052 0.05077 .
## c.mal -0.25958 0.28588 -0.908 0.37255
## ln.vol_cub 0.35518 0.11691 3.038 0.00551 **
## c.con_art5 -0.02596 0.01872 -1.387 0.17764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 0.9827792)
##
## Null deviance: 35.823 on 28 degrees of freedom
## Residual deviance: 24.749 on 25 degrees of freedom
## (44 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
res.artA <- resid(glm(rich_art_5 ~ ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1.2, subset = no_tim > 0 & host=="A" ), 'deviance')
res.malA <- resid(lm(c.mal ~ ln.vol_cub + c.con_art5,
data=n1.2, subset = no_tim > 0 & host=="A"), "working")
an.resA <- lm(res.artA ~ res.malA); summary(an.resA)
##
## Call:
## lm(formula = res.artA ~ res.malA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6110 -0.7195 -0.1231 0.5514 1.8781
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.09272 0.17442 -0.532 0.599
## res.malA -0.60915 0.46384 -1.313 0.200
##
## Residual standard error: 0.9393 on 27 degrees of freedom
## Multiple R-squared: 0.06004, Adjusted R-squared: 0.02523
## F-statistic: 1.725 on 1 and 27 DF, p-value: 0.2001
ndA=data.frame(res.malA=seq(min(res.malA), max(res.malA), len=50))
predsA <- predict(an.resA, newdata=ndA, type = "response", se.fit=TRUE)
ndA$fit <- predsA$fit
ndA$lower <- predsA$fit - predsA$se.fit
ndA$upper <- predsA$fit + predsA$se.fit
# Ceanothus
an3.fC <- glm(rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1, subset=host=="C")
summary(an3.fC)
##
## Call:
## glm(formula = rich_art_5 ~ c.mal + ln.vol_cub + c.con_art5, family = quasipoisson,
## data = n1, subset = host == "C")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3858 -1.1404 -0.2420 0.8255 2.0197
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.75955 0.79599 -3.467 0.00098 ***
## c.mal 0.72862 0.22775 3.199 0.00220 **
## ln.vol_cub 0.35643 0.06266 5.688 4.05e-07 ***
## c.con_art5 0.05688 0.01992 2.855 0.00590 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.538553)
##
## Null deviance: 184.732 on 63 degrees of freedom
## Residual deviance: 97.756 on 60 degrees of freedom
## (9 observations deleted due to missingness)
## AIC: NA
##
## Number of Fisher Scoring iterations: 5
res.artC <- resid(glm(rich_art_5 ~ ln.vol_cub + c.con_art5,
family=quasipoisson, data=n1.2, subset = no_tim > 0 & host=="C" ), 'deviance')
res.malC <- resid(lm(c.mal ~ ln.vol_cub + c.con_art5,
data=n1.2, subset = no_tim > 0 & host=="C"), "working")
an.resC <- lm(res.artC ~ res.malC); summary(an.resC)
##
## Call:
## lm(formula = res.artC ~ res.malC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.53433 -0.94271 -0.05687 1.10582 2.21744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1589 0.1580 -1.006 0.31827
## res.malC 1.3378 0.4632 2.888 0.00533 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.264 on 62 degrees of freedom
## Multiple R-squared: 0.1186, Adjusted R-squared: 0.1044
## F-statistic: 8.342 on 1 and 62 DF, p-value: 0.005327
ndC=data.frame(res.malC=seq(min(res.malC), max(res.malC), len=50))
predsC <- predict(an.resC, newdata=ndC, type = "response", se.fit=TRUE)
ndC$fit <- predsC$fit
ndC$lower <- predsC$fit - predsC$se.fit
ndC$upper <- predsC$fit + predsC$se.fit
jpeg(file="figures/n1_rich_mal.jpg",
width=750, height=750)
par(mar=c(6, 6, 1, 1))
plot(y=c(res.artA, res.artC), x=c(res.malA, res.malC),
xlab="", ylab="",
ylim=c(min(c(res.artA, res.artC)), 3),
las=1, type="n",
cex=2, cex.lab=2, cex.axis=2)
title(ylab="arthropod richness (residuals)", cex.lab=3, line=3.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4)
#axis.break(axis=2, breakpos=3.5)
points(y=res.artA, x=res.malA,
pch=16, col="blue", cex=2)
points(y=res.artC, x=res.malC,
pch=16, col="orange", cex=2)
polygon(x = c(ndC$res.malC, rev(ndC$res.malC)),
y = c(ndC$upper, rev(ndC$lower)),
col="orange", density=100, lwd=.5)
lines(ndC$res.malC, ndC$fit, lwd=2, col="orange")
legend(x=-.5, y=3,
legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))),
pch=c(16, 16), col=c("orange", "blue"),
bg=NULL, bty="n", cex=c(2, 2), y.intersp = 1)
dev.off()
## quartz_off_screen
## 2
n1.3 <- n1[complete.cases(n1$mas, n1$c.no_art_5, n1$c.mal), ]
# Adenostoma
an.s2A <- lm(mas ~ c.mal + c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
summary(an.s2A) # R-sq = 18.4 %
##
## Call:
## lm(formula = mas ~ c.mal + c.no_art_5, data = n1.3[n1.3$no_art_5 >
## 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62088 -0.20001 0.03397 0.15980 0.85752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.40958 0.10117 -4.048 0.000536 ***
## c.mal 0.18082 0.21821 0.829 0.416207
## c.no_art_5 -0.04825 0.01169 -4.127 0.000443 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3531 on 22 degrees of freedom
## Multiple R-squared: 0.4409, Adjusted R-squared: 0.3901
## F-statistic: 8.675 on 2 and 22 DF, p-value: 0.001668
an.mas.artA <- lm(mas ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
an.mal.artA <- lm(c.mal ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "A", ])
res.masA <- resid(an.mas.artA)
res.malA <- resid(an.mal.artA)
an.p_malA <- lm(res.masA ~ res.malA)
summary(an.p_malA) # r-sq = 13.5! on Adenostoma
##
## Call:
## lm(formula = res.masA ~ res.malA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62088 -0.20001 0.03397 0.15980 0.85752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.775e-18 6.906e-02 0.000 1.000
## res.malA 1.808e-01 2.134e-01 0.847 0.406
##
## Residual standard error: 0.3453 on 23 degrees of freedom
## Multiple R-squared: 0.03027, Adjusted R-squared: -0.0119
## F-statistic: 0.7179 on 1 and 23 DF, p-value: 0.4056
ndA=data.frame(res.malA=seq(min(res.malA), max(res.malA), len=50))
predsA <- predict(an.p_malA, newdata=ndA, type = "response", se.fit=TRUE)
ndA$fit <- predsA$fit
ndA$lower <- predsA$fit - predsA$se.fit
ndA$upper <- predsA$fit + predsA$se.fit
# Ceanothus
an.s2C <- lm(mas ~ c.mal + c.no_art_5 , data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
summary(an.s2C) # R-sq = 18.4 % (8.4)
##
## Call:
## lm(formula = mas ~ c.mal + c.no_art_5, data = n1.3[n1.3$no_art_5 >
## 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52193 -0.21290 -0.07037 0.18720 0.83286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.720632 0.047085 -15.305 < 2e-16 ***
## c.mal 0.113283 0.129225 0.877 0.385
## c.no_art_5 -0.029155 0.004516 -6.456 5.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2882 on 48 degrees of freedom
## Multiple R-squared: 0.465, Adjusted R-squared: 0.4427
## F-statistic: 20.86 on 2 and 48 DF, p-value: 3.025e-07
an.mas.artC <- lm(mas ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
an.mal.artC <- lm(c.mal ~ c.no_art_5, data=n1.3[n1.3$no_art_5 > 2 & n1.3$no_tim > 0 & n1.3$host == "C", ])
res.masC <- resid(an.mas.artC)
res.malC <- resid(an.mal.artC)
an.resC <- lm(res.masC ~ res.malC); summary(an.resC)
##
## Call:
## lm(formula = res.masC ~ res.malC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52193 -0.21290 -0.07037 0.18720 0.83286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.999e-17 3.994e-02 0.000 1.00
## res.malC 1.133e-01 1.279e-01 0.886 0.38
##
## Residual standard error: 0.2852 on 49 degrees of freedom
## Multiple R-squared: 0.01576, Adjusted R-squared: -0.004329
## F-statistic: 0.7845 on 1 and 49 DF, p-value: 0.3801
ndC=data.frame(res.malC=seq(min(res.malC), max(res.malC), len=50))
predsC <- predict(an.resC, newdata=ndC, type = "response", se.fit=TRUE)
ndC$fit <- predsC$fit
ndC$lower <- predsC$fit - predsC$se.fit
ndC$upper <- predsC$fit + predsC$se.fit
jpeg(file="figures/n1_mas_mal.jpg",
width=750, height=750)
par(mar=c(6, 7, 1, 1))
plot(y=c(res.masA, res.masC),
x=c(res.malA, res.malC),
xlab = "",
ylab = "",
las=1, type="n", cex.lab=2, cex.axis=2, cex=2)
title(ylab="mass-abundance slope (residuals)", cex.lab=3, line=4.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4)
points(y=res.masC, x=res.malC,
pch=16, col="orange", cex=2)
points(y=res.masA, x=res.malA,
pch=16, col="blue", cex=2)
#lines(nd$res.mal, nd$upper, lty=2)
#lines(nd$res.mal, nd$lower, lty=2)
polygon(x = c(ndA$res.malA, rev(ndA$res.malA)),
y = c(ndA$upper, rev(ndA$lower)),
col="blue", density=50, lwd=.5, border=NA)
lines(ndA$res.malA, ndA$fit,
col="blue", lwd=2)
legend(x=-.5, y=.35,
legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))),
pch=c(16, 16), col=c("orange", "blue"),
bg=NULL, bty="n", cex=2, y.intersp = 1)
dev.off()
## quartz_off_screen
## 2
p.n1 <- n1[complete.cases(n1$host, n1$c.mal, n1$cn2) & n1$pN < 1, ]
p.n1$c.host <- ifelse(p.n1$host == "A", -.5, .5)
an.mal <- lm(c.mal ~ c.host, data=p.n1)
an.cn_h <- lm(cn2 ~ c.host, data=p.n1)
r.mal <- resid(an.mal)
r.cn_h <- resid(an.cn_h)
an.res <- lm(r.cn_h ~ r.mal); summary(an.res)
##
## Call:
## lm(formula = r.cn_h ~ r.mal)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3758 -1.3005 -0.0045 1.1843 4.6236
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.601e-17 1.950e-01 0.000 1.0000
## r.mal -1.301e+00 5.417e-01 -2.403 0.0185 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.798 on 83 degrees of freedom
## Multiple R-squared: 0.06502, Adjusted R-squared: 0.05376
## F-statistic: 5.772 on 1 and 83 DF, p-value: 0.01851
nd=data.frame(r.mal=seq(min(r.mal), max(r.mal), len=50))
preds <- predict(an.res, newdata=nd, type = "response", se.fit=TRUE)
nd$fit <- preds$fit
nd$lower <- preds$fit - preds$se.fit
nd$upper <- preds$fit + preds$se.fit
jpeg(file="figures/n1_cn_mal.jpg",
width=750, height=750)
par(mar=c(6, 6, 1, 1))
plot(y=r.cn_h, x=r.mal,
xlab="", ylab="",
las=1, pch=16,
cex=2, cex.lab=3, cex.axis=2)
title(ylab="carbon:nitrogen ratio (residuals)", cex.lab=3, line=3.5)
title(xlab="maladaptation (residuals)", cex.lab=3, line=4)
#axis.break(axis=2, breakpos=3.5)
polygon(x = c(nd$r.mal, rev(nd$r.mal)),
y = c(nd$upper, rev(nd$lower)),
col="grey", density=100, lwd=.5)
lines(nd$r.mal, nd$fit, lwd=2)
dev.off()
## quartz_off_screen
## 2
# carbon on CN
jpeg(file="figures/n1_cn_pC.jpg",
width=750, height=750)
par(mar=c(6, 7, 3, 1))
plot(y=p.n1$cn2, x=p.n1$pC,
ylab="", xlab="",
las=1, pch=16, col=ifelse(p.n1$host == "A", "blue", "orange"),
cex=2, cex.lab=3, cex.axis=2)
title(ylab="carbon:nitrogen ratio", line=4.5, cex.lab=3)
title(xlab="mass carbon (mg)", line=4.5, cex.lab=3)
legend(x=7.25, y=33,
legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))),
pch=c(16, 16), col=c("orange", "blue"),
bg=NULL, bty="n", cex=2, y.intersp = 1)
dev.off()
## quartz_off_screen
## 2
# nitrogen on CN
an_n <- lm(cn2 ~ pN + I(pN ^ 2), data=p.n1)
summary(an_n)
##
## Call:
## lm(formula = cn2 ~ pN + I(pN^2), data = p.n1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1372 -0.9234 0.1656 1.0417 4.3042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.181 5.368 11.584 < 2e-16 ***
## pN -133.353 24.095 -5.535 3.65e-07 ***
## I(pN^2) 100.705 26.418 3.812 0.000266 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.828 on 82 degrees of freedom
## Multiple R-squared: 0.8013, Adjusted R-squared: 0.7964
## F-statistic: 165.3 on 2 and 82 DF, p-value: < 2.2e-16
nd <- data.frame(pN=seq(min(p.n1$pN), max(p.n1$pN), len=50))
preds <- predict(an_n, newdata=nd, type="response", se.fit=TRUE)
nd <- cbind(nd, fit=preds$fit, upper=preds$fit + preds$se.fit, lower=preds$fit-preds$se.fit)
jpeg(file="figures/n1_cn_pN.jpg",
width=750, height=750)
par(mar=c(6, 7, 1, 1))
plot(y=p.n1$cn2, x=p.n1$pN,
ylab="", xlab="",
las=1, pch=16, col=ifelse(p.n1$host == "A", "blue", "orange"),
cex=2, cex.lab=2, cex.axis=2)
title(ylab="carbon:nitrogen ratio", line=4.5, cex.lab=3)
title(xlab="mass nitrogen (mg)", line=4.5, cex.lab=3)
polygon(x = c(nd$pN, rev(nd$pN)),
y = c(nd$upper, rev(nd$lower)),
col="grey", density=100, lwd=.5)
lines(nd$pN, nd$fit, lwd=2, col="grey")
legend(x=.5, y=30,
legend=c(expression(italic("C. spinosus"), italic( "A. fasciculatum"))),
pch=c(16, 16), col=c("orange", "blue"),
bg=NULL, bty="n", cex=2, y.intersp = 1)
dev.off()
## quartz_off_screen
## 2
library(scales)
## Warning: package 'scales' was built under R version 3.6.2
##
## Attaching package: 'scales'
## The following object is masked from 'package:plotrix':
##
## rescale
n1$occ <- ifelse(n1$no_tim_gs==0, FALSE, TRUE)
n1$radius <- (n1$vol_cub/(100000*pi))^(1/3)
n1_occ <- n1[n1$occ,]
adOnly <- n1_occ[n1_occ$host=="A",]
ceOnly <- n1_occ[n1_occ$host=="C",]
empty <- n1[!n1$occ,]
cFunA <- colorRamp(c("skyblue1", "blue"))
cFunC <- colorRamp(c("darkorange2", "peachpuff"))
colsA <- rgb(cFunA(adOnly$mal), maxColorValue=255)
colsC <- rgb(cFunC(ceOnly$mal), maxColorValue=255)
jpeg(file="figures/n1_map.jpg",
width=750, height=750)
par(mar=c(1, 1, 1, 1))
plot(x=1:10, y=1:10,
xlim=c(0,80), ylim=c(-1, 50),
type="n",
xlab="", ylab="",
las=1, cex.axis=2, cex.lab=2, axes=FALSE)
#title(ylab="meters N-S", line=4, cex.lab=2)
polygon(x=c(-3, -3, 85, 85), y=c(-3, 55, 55, -3), density=500, col="lightgrey")
box(which="plot", lty="solid")
polygon(x=c(60, 60, 80, 80), y=c(0, .25, .25, 0), density=500, col="black")
text(70, -.75, labels="10 meters", cex=1.5)
symbols(x=adOnly$x, y=adOnly$y, circles=(adOnly$radius), bg=colsA, inches=FALSE,
xlim=c(0,75), ylim=c(-1, 50),
xlab = "metres E-W", ylab="metres N-S", fg="blue", lwd=2.5, cex.lab=1.5,
cex.axis=1.5, add=TRUE)
symbols(x=ceOnly$x, y=ceOnly$y, circles=(ceOnly$radius), bg=colsC, inches=F, add=T,
fg="orange", lwd=2.5)
symbols(x=empty$x, y=empty$y, circles=(empty$radius), bg="white",
fg=alpha(ifelse(empty$host=="A", "blue", "orange"),0.5),
inches=F, add=T, lwd=2.5)
legend(x=0, y=50, legend=c(expression(italic("A. fasciculatum")), expression(italic("C. spinosus"))), pch=16, col=c("blue", "orange"), bty="n", cex=2)
dev.off()
## quartz_off_screen
## 2
ad$taxon <- factor(ifelse(ad$family == "Formicidae", "Formicidae",
ifelse(ad$order == "Homoptera", "Hemiptera",
ifelse(ad$order == "Neuroptera", "Chrysopidae",
ifelse(ad$order == "Aranea",
as.character(ad$family),
ifelse(ad$order == "Coleoptera" & ad$family != "",
as.character(ad$family),
as.character(ad$order)))))))
#ad$taxon <- factor(ifelse(ad$family != "", as.character(ad$family),
# as.character(ad$order)))
ad$taxon <- factor(ifelse(ad$family == "Formicidae", "Formicidae",
ifelse(ad$order == "Homoptera", "Hemiptera",
ifelse(ad$order == "Aranea", "Aranea",
ifelse(ad$order %in% c("Chrysopidae", "Diptera",
"Hymenoptera"),
"Other", as.character(ad$order))))))
orders <- levels(droplevels(ad$taxon[ad$length >=5]))
arts <- setNames(lapply(orders,
FUN=function(x) table(ad$plant_id[ad$taxon == x & ad$length >=5])),
orders)
newnames <- paste("no_", orders, "_5", sep='')
n1[, newnames] <- 0
for(i in 1:length(newnames)) {
n1[match(names(arts[[i]]), n1$plant_id), newnames[i]] <- arts[[i]]
}
pnames <- paste("p_", orders, sep="")
n1[, pnames] <- 0
for(i in 1:length(pnames)) {
n1[, pnames[i]] <- ifelse(n1$no_any_5 == 0, NA, n1[, newnames[i]] / n1[, "no_any_5"])
}
thresh <- c(seq(0, .9, length.out = 10), 0.9999999)
surv <- do.call(rbind, lapply(pnames, FUN = function(x) {
sapply(thresh, FUN=function(y) {
sum(n1[, x] > y, na.rm=TRUE)
})
})
)
dimnames(surv) <- list(orders, c(thresh[1:10], 1))
arttot <- table(droplevels(ad$taxon[ad$length >= 5]))
surv <- cbind(tot=arttot, surv)
surv <- surv[order(surv[, "tot"], decreasing=TRUE), ]
#ad[ad$order == "Lepidoptera" & ad$length >= 5, c("order", "family")]
#write.csv(surv, file="~/Dropbox/Projects/CA Arthropods/N1/Data/art_surv.csv",
# row.names=TRUE)
write.csv(surv, file="data/art_surv_gp.csv",
row.names=TRUE)
# smaller things
arttot <- table(droplevels(ad$taxon[ad$length <5]))
ad[ad$family=="Formicidae", c("genus")]
## NULL
library(partitions)
p = 0.05
s = 10
naive.pmf <- function(prob = .05, giveup = 10, max = 40, plot.pmf = TRUE) {
p = prob
s = giveup
pmf <- c((1 - pgeom(s-1, p)), c(lapply(1:(max-s), function(x) {
library(partitions)
part <- t(parts(x))
ps <- lapply(
lapply(1:nrow(part), function(y) part[y,]), function(z) z[z > 0])
ps <- ps[sapply(ps, function(q) all(q < s))]
pt.ct <- mapply(FUN=function(unq, len) factorial(len)/prod(factorial(unq)),
unq=lapply(ps, table),
len=lapply(ps, length))
probs <- mapply(function(ps, pt.ct) prod(dgeom(ps-1, prob=p)) * pt.ct,
ps = ps, pt.ct = pt.ct)
sum(probs) * (1 - pgeom(s-1, p))
}), recursive=TRUE))
if(plot.pmf == TRUE) plot(10:max, pmf, type="l")
cum <- sum(pmf) # cumulative density = 0.948273
expt <- sum(pmf * 10:max) # expected value = 12.0177
list(pmf = pmf, cum = cum, expt = expt)
}
dgeom(0, .05)
## [1] 0.05
dgeom(0, .05)^2 + dgeom(1, .05)
## [1] 0.05
dgeom(0, .05)^3 + 2*(dgeom(1, .05) * dgeom(0, .05)) + dgeom(2, .05)
## [1] 0.05
lm(c(.3, .65) ~ c(.2, .8))
##
## Call:
## lm(formula = c(0.3, 0.65) ~ c(0.2, 0.8))
##
## Coefficients:
## (Intercept) c(0.2, 0.8)
## 0.1833 0.5833
st.w <- function(freq) { 1.117 - 0.583 * freq }
gr.w <- function(freq) { 0.183 + 0.583 * freq}
w.bar <- function(freq) st.w(freq) * freq + gr.w(freq) * (1 - freq)
wbs <- w.bar(seq(0, 1, .01))
cbind(seq(0,1,.01), wbs)
## wbs
## [1,] 0.00 0.1830000
## [2,] 0.01 0.1980534
## [3,] 0.02 0.2128736
## [4,] 0.03 0.2274606
## [5,] 0.04 0.2418144
## [6,] 0.05 0.2559350
## [7,] 0.06 0.2698224
## [8,] 0.07 0.2834766
## [9,] 0.08 0.2968976
## [10,] 0.09 0.3100854
## [11,] 0.10 0.3230400
## [12,] 0.11 0.3357614
## [13,] 0.12 0.3482496
## [14,] 0.13 0.3605046
## [15,] 0.14 0.3725264
## [16,] 0.15 0.3843150
## [17,] 0.16 0.3958704
## [18,] 0.17 0.4071926
## [19,] 0.18 0.4182816
## [20,] 0.19 0.4291374
## [21,] 0.20 0.4397600
## [22,] 0.21 0.4501494
## [23,] 0.22 0.4603056
## [24,] 0.23 0.4702286
## [25,] 0.24 0.4799184
## [26,] 0.25 0.4893750
## [27,] 0.26 0.4985984
## [28,] 0.27 0.5075886
## [29,] 0.28 0.5163456
## [30,] 0.29 0.5248694
## [31,] 0.30 0.5331600
## [32,] 0.31 0.5412174
## [33,] 0.32 0.5490416
## [34,] 0.33 0.5566326
## [35,] 0.34 0.5639904
## [36,] 0.35 0.5711150
## [37,] 0.36 0.5780064
## [38,] 0.37 0.5846646
## [39,] 0.38 0.5910896
## [40,] 0.39 0.5972814
## [41,] 0.40 0.6032400
## [42,] 0.41 0.6089654
## [43,] 0.42 0.6144576
## [44,] 0.43 0.6197166
## [45,] 0.44 0.6247424
## [46,] 0.45 0.6295350
## [47,] 0.46 0.6340944
## [48,] 0.47 0.6384206
## [49,] 0.48 0.6425136
## [50,] 0.49 0.6463734
## [51,] 0.50 0.6500000
## [52,] 0.51 0.6533934
## [53,] 0.52 0.6565536
## [54,] 0.53 0.6594806
## [55,] 0.54 0.6621744
## [56,] 0.55 0.6646350
## [57,] 0.56 0.6668624
## [58,] 0.57 0.6688566
## [59,] 0.58 0.6706176
## [60,] 0.59 0.6721454
## [61,] 0.60 0.6734400
## [62,] 0.61 0.6745014
## [63,] 0.62 0.6753296
## [64,] 0.63 0.6759246
## [65,] 0.64 0.6762864
## [66,] 0.65 0.6764150
## [67,] 0.66 0.6763104
## [68,] 0.67 0.6759726
## [69,] 0.68 0.6754016
## [70,] 0.69 0.6745974
## [71,] 0.70 0.6735600
## [72,] 0.71 0.6722894
## [73,] 0.72 0.6707856
## [74,] 0.73 0.6690486
## [75,] 0.74 0.6670784
## [76,] 0.75 0.6648750
## [77,] 0.76 0.6624384
## [78,] 0.77 0.6597686
## [79,] 0.78 0.6568656
## [80,] 0.79 0.6537294
## [81,] 0.80 0.6503600
## [82,] 0.81 0.6467574
## [83,] 0.82 0.6429216
## [84,] 0.83 0.6388526
## [85,] 0.84 0.6345504
## [86,] 0.85 0.6300150
## [87,] 0.86 0.6252464
## [88,] 0.87 0.6202446
## [89,] 0.88 0.6150096
## [90,] 0.89 0.6095414
## [91,] 0.90 0.6038400
## [92,] 0.91 0.5979054
## [93,] 0.92 0.5917376
## [94,] 0.93 0.5853366
## [95,] 0.94 0.5787024
## [96,] 0.95 0.5718350
## [97,] 0.96 0.5647344
## [98,] 0.97 0.5574006
## [99,] 0.98 0.5498336
## [100,] 0.99 0.5420334
## [101,] 1.00 0.5340000
max(wbs)
## [1] 0.676415
plot(seq(0, 1, .01), wbs)
#wss <- ws(seq(0, 1, .001))
#max(wss - min(wss))
#min(wss)
ws <- function(freq) {
raw <- 1.517 * freq - 1.166 * freq^2 + 0.183
(raw - 0.183)/.493415
}
ws <- function(freq) {
raw <- w.bar(freq)
(raw - 0.183)/.493415
}
ws(seq(0, 1, .01))
## [1] 0.00000000 0.03050860 0.06054457 0.09010792 0.11919865 0.14781675
## [7] 0.17596222 0.20363507 0.23083530 0.25756290 0.28381788 0.30960024
## [13] 0.33490996 0.35974707 0.38411155 0.40800340 0.43142264 0.45436924
## [19] 0.47684323 0.49884458 0.52037332 0.54142943 0.56201291 0.58212377
## [25] 0.60176201 0.62092762 0.63962060 0.65784097 0.67558870 0.69286382
## [31] 0.70966631 0.72599617 0.74185341 0.75723802 0.77215002 0.78658938
## [37] 0.80055612 0.81405024 0.82707173 0.83962060 0.85169685 0.86330047
## [43] 0.87443146 0.88508983 0.89527558 0.90498870 0.91422920 0.92299707
## [49] 0.93129232 0.93911494 0.94646494 0.95334232 0.95974707 0.96567919
## [55] 0.97113870 0.97612557 0.98063983 0.98468145 0.98825046 0.99134684
## [61] 0.99397059 0.99612172 0.99780023 0.99900611 0.99973937 1.00000000
## [67] 0.99978801 0.99910339 0.99794615 0.99631629 0.99421380 0.99163868
## [73] 0.98859094 0.98507058 0.98107759 0.97661198 0.97167374 0.96626288
## [79] 0.96037940 0.95402329 0.94719455 0.93989319 0.93211921 0.92387260
## [85] 0.91515337 0.90596151 0.89629703 0.88615993 0.87555020 0.86446784
## [91] 0.85291286 0.84088526 0.82838503 0.81541218 0.80196670 0.78804860
## [97] 0.77365787 0.75879452 0.74345855 0.72764995 0.71136873
plot(seq(0, 1, .01), ws(seq(0, 1, .01)))
freq <- seq(0, .99, .01)
k = 1
ws <- 1- (1 / (1 - exp(freq)))
plot(freq, ws)
freq <- c(0.2, 0.8)
df <- data.frame(a=c(.99, .65), b=I(1 - exp(-freq)))
an1 <- lm(a ~ b, data=df)
preds <- predict(object = an1, newdata = data.frame(b=seq(0, .99, .01)))
k = 1
c = .9
min.RSS <- function(data, par) {
with(data, sum((par[1] + par[2] * (1 / (1 - c * exp(-x))) - y)^2))
}
dat <- data.frame(x = c(0.2, 0.8), y = c(1, .65))
pars <- optim(par = c(0, 1, 1), min.RSS, data=dat)$par
freqs <- seq(0, .99, .01)
ws <- pars[1] + pars[2] * (1 / (1 - c * exp(-freqs)))
plot(freqs, ws)
nfmal <- function(mal=NULL, k=1) {
freq <- 1 - mal
# adapted moph
min.RSS <- function(data, par) {
with(data, sum(( (par[1] / (1 - ((par[2] - par[1])/par[2]) * exp((-1) * k * x))) - y)^2))
}
dat.ad <- data.frame(x = c(0.2, 0.8), y = c(1, .65))
pars.ad <- optim(par = c(0,1), min.RSS, data=dat.ad)$par
freqs <- seq(0, 1, .01)
ws.ad <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freqs)))
# plot(freqs, ws.ad)
# maladapted morph
ws.mal <- rev(ws.ad) - .35
raw <- (ws.ad * freqs) + (ws.mal * (1 - freqs))
min.w <- min(raw)
max.w <- max(raw - min.w)
scaled <- (raw - min.w)/max.w
cbind(freqs, scaled)
plot(rev(freqs), 1 - scaled)
#plot(freqs, ws.ad, ylim=c(-2, 2), type="l")
#plot(freqs, ws.ad, type="l")
#lines(freqs, ws.mal, lty=2)
# for output
ws.ad2 <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freq)))
ws.mal2 <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * (1-freq)))) - .35
raw2 <- (ws.ad2 * freq) + (ws.mal2 * (1 - freq))
sc2 <- (raw2-min.w) / max.w
#plot(freqs, mal)
#pars.ad
1 - sc2
}
#lm(c(1, .65) ~ c(.16, ))
st.w <- function(freq) { 1.117 - 0.583 * freq*(1-freq)}
gr.w <- function(freq) { 0.183 + 0.583 * freq*(1-freq)}
w.bar <- function(freq) st.w(freq-.01) * freq + gr.w(freq-.01) * (1 - (freq-.01))
wbs <- w.bar(seq(0, 1, .01))
cbind(seq(0,1,.01), wbs)
## wbs
## [1,] 0.00 0.1788828
## [2,] 0.01 0.1941700
## [3,] 0.02 0.2091085
## [4,] 0.03 0.2237055
## [5,] 0.04 0.2379677
## [6,] 0.05 0.2519024
## [7,] 0.06 0.2655163
## [8,] 0.07 0.2788166
## [9,] 0.08 0.2918103
## [10,] 0.09 0.3045043
## [11,] 0.10 0.3169056
## [12,] 0.11 0.3290213
## [13,] 0.12 0.3408583
## [14,] 0.13 0.3524236
## [15,] 0.14 0.3637242
## [16,] 0.15 0.3747672
## [17,] 0.16 0.3855594
## [18,] 0.17 0.3961080
## [19,] 0.18 0.4064198
## [20,] 0.19 0.4165020
## [21,] 0.20 0.4263615
## [22,] 0.21 0.4360052
## [23,] 0.22 0.4454402
## [24,] 0.23 0.4546735
## [25,] 0.24 0.4637121
## [26,] 0.25 0.4725630
## [27,] 0.26 0.4812331
## [28,] 0.27 0.4897295
## [29,] 0.28 0.4980592
## [30,] 0.29 0.5062291
## [31,] 0.30 0.5142463
## [32,] 0.31 0.5221177
## [33,] 0.32 0.5298504
## [34,] 0.33 0.5374513
## [35,] 0.34 0.5449274
## [36,] 0.35 0.5522858
## [37,] 0.36 0.5595334
## [38,] 0.37 0.5666773
## [39,] 0.38 0.5737243
## [40,] 0.39 0.5806816
## [41,] 0.40 0.5875561
## [42,] 0.41 0.5943548
## [43,] 0.42 0.6010847
## [44,] 0.43 0.6077528
## [45,] 0.44 0.6143661
## [46,] 0.45 0.6209316
## [47,] 0.46 0.6274563
## [48,] 0.47 0.6339472
## [49,] 0.48 0.6404113
## [50,] 0.49 0.6468555
## [51,] 0.50 0.6532869
## [52,] 0.51 0.6597125
## [53,] 0.52 0.6661392
## [54,] 0.53 0.6725742
## [55,] 0.54 0.6790242
## [56,] 0.55 0.6854965
## [57,] 0.56 0.6919978
## [58,] 0.57 0.6985353
## [59,] 0.58 0.7051160
## [60,] 0.59 0.7117468
## [61,] 0.60 0.7184347
## [62,] 0.61 0.7251868
## [63,] 0.62 0.7320100
## [64,] 0.63 0.7389113
## [65,] 0.64 0.7458977
## [66,] 0.65 0.7529763
## [67,] 0.66 0.7601539
## [68,] 0.67 0.7674377
## [69,] 0.68 0.7748345
## [70,] 0.69 0.7823515
## [71,] 0.70 0.7899956
## [72,] 0.71 0.7977737
## [73,] 0.72 0.8056929
## [74,] 0.73 0.8137602
## [75,] 0.74 0.8219826
## [76,] 0.75 0.8303671
## [77,] 0.76 0.8389206
## [78,] 0.77 0.8476502
## [79,] 0.78 0.8565629
## [80,] 0.79 0.8656656
## [81,] 0.80 0.8749654
## [82,] 0.81 0.8844692
## [83,] 0.82 0.8941841
## [84,] 0.83 0.9041170
## [85,] 0.84 0.9142749
## [86,] 0.85 0.9246649
## [87,] 0.86 0.9352939
## [88,] 0.87 0.9461690
## [89,] 0.88 0.9572970
## [90,] 0.89 0.9686851
## [91,] 0.90 0.9803402
## [92,] 0.91 0.9922693
## [93,] 0.92 1.0044794
## [94,] 0.93 1.0169775
## [95,] 0.94 1.0297706
## [96,] 0.95 1.0428657
## [97,] 0.96 1.0562698
## [98,] 0.97 1.0699899
## [99,] 0.98 1.0840330
## [100,] 0.99 1.0984060
## [101,] 1.00 1.1131160
max(wbs)
## [1] 1.113116
plot(seq(0, 1, .01), wbs)
#wss <- ws(seq(0, 1, .001))
#max(wss - min(wss))
#min(wss)
ws <- function(freq) {
w.bar <- function(freq) st.w(freq-.01) * freq + gr.w(freq-.01) * (1 - freq-.01)
raw <- w.bar(freq) * freq
raw/1.10934
}
ws(seq(0, 1, .01))
## [1] 0.000000000 0.001717327 0.003701897 0.005944533 0.008436311 0.011168560
## [7] 0.014132861 0.017321048 0.020725205 0.024337671 0.028151034 0.032158138
## [13] 0.036352077 0.040726197 0.045274098 0.049989630 0.054866897 0.059900256
## [19] 0.065084313 0.070413929 0.075884216 0.081490540 0.087228518 0.093094017
## [25] 0.099083161 0.105192323 0.111418129 0.117757457 0.124207438 0.130765455
## [31] 0.137429142 0.144196388 0.151065331 0.158034364 0.165102130 0.172267526
## [37] 0.179529701 0.186888056 0.194342243 0.201892168 0.209537989 0.217280116
## [43] 0.225119211 0.233056187 0.241092213 0.249228707 0.257467340 0.265810035
## [49] 0.274258969 0.282816569 0.291485515 0.300268741 0.309169430 0.318191021
## [55] 0.327337201 0.336611913 0.346019350 0.355563959 0.365250438 0.375083737
## [61] 0.385069059 0.395211859 0.405517845 0.415992976 0.426643465 0.437475774
## [67] 0.448496621 0.459712974 0.471132054 0.482761334 0.494608540 0.506681649
## [73] 0.518988892 0.531538749 0.544339957 0.557401501 0.570732620 0.584342806
## [79] 0.598241803 0.612439605 0.626946460 0.641772870 0.656929587 0.672427615
## [85] 0.688278211 0.704492885 0.721083398 0.738061764 0.755440249 0.773231371
## [91] 0.791447901 0.810102862 0.829209529 0.848781429 0.868832341 0.889376298
## [97] 0.910427584 0.932000735 0.954110539 0.976772039 1.000000526
plot(seq(0, 1, .01), ws(seq(0, 1, .01)))
freqs <- seq(0, 1, .01)
png(filename = "figures/nfds_mal.png")
plot(freqs, freqs, type="n",
xlab= "frequency of cryptic morph",
ylab= "average population fitness (scaled)",
cex.lab = 1.5, las = 1, cex.axis = 1.5)
for (k in c(0.01, .6 , 1.2, 2, 3, 4)) {
min.RSS <- function(data, par) {
with(data, sum(( (par[1] / (1 - ((par[2] - par[1])/par[2]) * exp((-1) * k * x))) - y)^2))
}
dat.ad <- data.frame(x = c(0.2, 0.8), y = c(1, .65))
pars.ad <- optim(par = c(0,1), min.RSS, data=dat.ad)$par
freqs <- seq(0, 1, .01)
ws.ad <- (pars.ad[1] / (1 - ((pars.ad[2] - pars.ad[1])/pars.ad[2]) * exp((-1) * k * freqs)))
# maladapted morph
ws.mal <- rev(ws.ad) - .35
raw <- (ws.ad * freqs) + (ws.mal * (1 - freqs))
min.w <- min(raw)
max.w <- max(raw - min.w)
scaled <- (raw - min.w)/max.w
lines(freqs, scaled)
}
text(c(.2, .7), c(.8, .6),
labels=c("k = 0.01", "k = 4.0"),
cex=1.5)
dev.off()
## quartz_off_screen
## 2
lm(c(1.667, 1.0835) ~ c(.2, .8))
##
## Call:
## lm(formula = c(1.667, 1.0835) ~ c(0.2, 0.8))
##
## Coefficients:
## (Intercept) c(0.2, 0.8)
## 1.8615 -0.9725
lm(c(.5, 1.0835) ~ c(.2, .8))
##
## Call:
## lm(formula = c(0.5, 1.0835) ~ c(0.2, 0.8))
##
## Coefficients:
## (Intercept) c(0.2, 0.8)
## 0.3055 0.9725
st.w <- function(freq) { 1.861 - .972 * (freq/(1-freq)) }
gr.w <- function(freq) { 0.306 + .972 * (freq/(1-freq)) }
w.bar <- function(freq) st.w(freq) * freq + gr.w(freq) * (1 - freq)
wbs <- w.bar(seq(0, .99, .01))
wbs.sc <- (wbs - min(wbs))/max(wbs-min(wbs))
cbind(seq(0,.99,.01), wbs.sc)
## wbs.sc
## [1,] 0.00 0.9916345
## [2,] 0.01 0.9919035
## [3,] 0.02 0.9921705
## [4,] 0.03 0.9924352
## [5,] 0.04 0.9926977
## [6,] 0.05 0.9929578
## [7,] 0.06 0.9932155
## [8,] 0.07 0.9934706
## [9,] 0.08 0.9937232
## [10,] 0.09 0.9939732
## [11,] 0.10 0.9942203
## [12,] 0.11 0.9944647
## [13,] 0.12 0.9947060
## [14,] 0.13 0.9949443
## [15,] 0.14 0.9951795
## [16,] 0.15 0.9954114
## [17,] 0.16 0.9956399
## [18,] 0.17 0.9958649
## [19,] 0.18 0.9960863
## [20,] 0.19 0.9963039
## [21,] 0.20 0.9965176
## [22,] 0.21 0.9967272
## [23,] 0.22 0.9969326
## [24,] 0.23 0.9971337
## [25,] 0.24 0.9973302
## [26,] 0.25 0.9975219
## [27,] 0.26 0.9977087
## [28,] 0.27 0.9978904
## [29,] 0.28 0.9980668
## [30,] 0.29 0.9982375
## [31,] 0.30 0.9984025
## [32,] 0.31 0.9985614
## [33,] 0.32 0.9987140
## [34,] 0.33 0.9988600
## [35,] 0.34 0.9989991
## [36,] 0.35 0.9991309
## [37,] 0.36 0.9992552
## [38,] 0.37 0.9993715
## [39,] 0.38 0.9994795
## [40,] 0.39 0.9995788
## [41,] 0.40 0.9996690
## [42,] 0.41 0.9997495
## [43,] 0.42 0.9998199
## [44,] 0.43 0.9998796
## [45,] 0.44 0.9999281
## [46,] 0.45 0.9999648
## [47,] 0.46 0.9999890
## [48,] 0.47 1.0000000
## [49,] 0.48 0.9999970
## [50,] 0.49 0.9999793
## [51,] 0.50 0.9999458
## [52,] 0.51 0.9998958
## [53,] 0.52 0.9998280
## [54,] 0.53 0.9997415
## [55,] 0.54 0.9996349
## [56,] 0.55 0.9995070
## [57,] 0.56 0.9993563
## [58,] 0.57 0.9991811
## [59,] 0.58 0.9989798
## [60,] 0.59 0.9987505
## [61,] 0.60 0.9984909
## [62,] 0.61 0.9981989
## [63,] 0.62 0.9978718
## [64,] 0.63 0.9975069
## [65,] 0.64 0.9971008
## [66,] 0.65 0.9966502
## [67,] 0.66 0.9961511
## [68,] 0.67 0.9955991
## [69,] 0.68 0.9949892
## [70,] 0.69 0.9943158
## [71,] 0.70 0.9935725
## [72,] 0.71 0.9927523
## [73,] 0.72 0.9918467
## [74,] 0.73 0.9908463
## [75,] 0.74 0.9897402
## [76,] 0.75 0.9885157
## [77,] 0.76 0.9871580
## [78,] 0.77 0.9856497
## [79,] 0.78 0.9839702
## [80,] 0.79 0.9820952
## [81,] 0.80 0.9799953
## [82,] 0.81 0.9776350
## [83,] 0.82 0.9749709
## [84,] 0.83 0.9719493
## [85,] 0.84 0.9685033
## [86,] 0.85 0.9645479
## [87,] 0.86 0.9599741
## [88,] 0.87 0.9546390
## [89,] 0.88 0.9483525
## [90,] 0.89 0.9408549
## [91,] 0.90 0.9317830
## [92,] 0.91 0.9206120
## [93,] 0.92 0.9065547
## [94,] 0.93 0.8883742
## [95,] 0.94 0.8640088
## [96,] 0.95 0.8297477
## [97,] 0.96 0.7781690
## [98,] 0.97 0.6919551
## [99,] 0.98 0.5191534
## [100,] 0.99 0.0000000
max(wbs)
## [1] 1.088568
plot(seq(0, .99, .01), wbs.sc)
png(filename = "figures/NFDS_model.png")
#par(mfrow=c(2,1))
plot(x=c(0,1), y=c(0,1.2), type="n",
xlab="cryptic morph frequency",
ylab="absolute fitness",
las=1)
abline(a=1.117, b = -0.583)
abline(a=0.183, b = 0.583, lty=2 )
abline(h=.65, lty=3, lwd=.5)
legend("topright",
legend=c("cryptic", "conspicuous"),
lty=c(1,2),
bty="n", cex = .75)
#plot(seq(0, 1, .01), wbs.sc, type="l", xlab="cryptic morph frequency",ylab="scaled average absolute fitness", las=1)
dev.off()
## quartz_off_screen
## 2
